Merge branch 'master' of http://git.icms.temple.edu/lammps-ro
Conflicts: doc/Section_intro.txt doc/Section_start.txt src/GPU/Install.sh src/GPU/fix_gpu.cpp
This commit is contained in:
Mike Brown
@ -399,12 +399,13 @@ potentials. Click on the style itself for a full description:
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long/gpu</A></TD><TD ><A HREF = "pair_charmm.html">lj/charmm/coul/long/opt</A></TD><TD ><A HREF = "pair_class2.html">lj/class2</A></TD><TD ><A HREF = "pair_class2.html">lj/class2/coul/cut</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_class2.html">lj/class2/coul/long</A></TD><TD ><A HREF = "pair_lj.html">lj/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/gpu</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/opt</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/cut/gpu</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/debye</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/gpu</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/tip4p</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs</A></TD><TD ><A HREF = "pair_lj_smooth.html">lj/smooth</A></TD><TD ><A HREF = "pair_lj96_cut.html">lj96/cut</A></TD><TD ><A HREF = "pair_lj96_cut.html">lj96/cut/gpu</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_lubricate.html">lubricate</A></TD><TD ><A HREF = "pair_meam.html">meam</A></TD><TD ><A HREF = "pair_morse.html">morse</A></TD><TD ><A HREF = "pair_morse.html">morse/opt</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_peri.html">peri/lps</A></TD><TD ><A HREF = "pair_peri.html">peri/pmb</A></TD><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_resquared.html">resquared</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_soft.html">soft</A></TD><TD ><A HREF = "pair_sw.html">sw</A></TD><TD ><A HREF = "pair_table.html">table</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa</A></TD><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid</A>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/gpu</A></TD><TD ><A HREF = "pair_lj.html">lj/cut/coul/long/tip4p</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand</A></TD><TD ><A HREF = "pair_lj_expand.html">lj/expand/gpu</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_gromacs.html">lj/gromacs</A></TD><TD ><A HREF = "pair_gromacs.html">lj/gromacs/coul/gromacs</A></TD><TD ><A HREF = "pair_lj_smooth.html">lj/smooth</A></TD><TD ><A HREF = "pair_lj96_cut.html">lj96/cut</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_lj96_cut.html">lj96/cut/gpu</A></TD><TD ><A HREF = "pair_lubricate.html">lubricate</A></TD><TD ><A HREF = "pair_meam.html">meam</A></TD><TD ><A HREF = "pair_morse.html">morse</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_morse.html">morse/gpu</A></TD><TD ><A HREF = "pair_morse.html">morse/opt</A></TD><TD ><A HREF = "pair_peri.html">peri/lps</A></TD><TD ><A HREF = "pair_peri.html">peri/pmb</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_reax.html">reax</A></TD><TD ><A HREF = "pair_resquared.html">resquared</A></TD><TD ><A HREF = "pair_soft.html">soft</A></TD><TD ><A HREF = "pair_sw.html">sw</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_table.html">table</A></TD><TD ><A HREF = "pair_tersoff.html">tersoff</A></TD><TD ><A HREF = "pair_tersoff_zbl.html">tersoff/zbl</A></TD><TD ><A HREF = "pair_yukawa.html">yukawa</A></TD></TR>
|
||||
<TR ALIGN="center"><TD ><A HREF = "pair_yukawa_colloid.html">yukawa/colloid</A>
|
||||
</TD></TR></TABLE></DIV>
|
||||
|
||||
<P>These are pair styles contributed by users, which can be used if
|
||||
@ -483,7 +484,8 @@ description:
|
||||
Kspace solvers. Click on the style itself for a full description:
|
||||
</P>
|
||||
<DIV ALIGN=center><TABLE BORDER=1 >
|
||||
<TR ALIGN="center"><TD WIDTH="100"><A HREF = "kspace_style.html">ewald</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm/tip4p</A>
|
||||
<TR ALIGN="center"><TD WIDTH="100"><A HREF = "kspace_style.html">ewald</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm/gpu/single</A></TD><TD WIDTH="100"><A HREF = "kspace_style.html">pppm/gpu/double</A></TD></TR>
|
||||
<TR ALIGN="center"><TD WIDTH="100"><A HREF = "kspace_style.html">pppm/tip4p</A>
|
||||
</TD></TR></TABLE></DIV>
|
||||
|
||||
<P>These are Kspace solvers contributed by users, which can be used if
|
||||
|
||||
@ -173,6 +173,10 @@ the bond topologies you have defined.
|
||||
neighbors for each atom. This likely means something is wrong with
|
||||
the bond topologies you have defined.
|
||||
|
||||
<DT><I>Accelerated style in input script but no fix gpu</I>
|
||||
|
||||
<DD>GPU acceleration requires fix gpu in the input script.
|
||||
|
||||
<DT><I>All angle coeffs are not set</I>
|
||||
|
||||
<DD>All angle coefficients must be set in the data file or by the
|
||||
@ -1240,9 +1244,9 @@ non-periodic z dimension.
|
||||
unless you use the kspace_modify command to define a 2d slab with a
|
||||
non-periodic z dimension.
|
||||
|
||||
<DT><I>Cannot use pair hybrid with multiple GPU pair styles</I>
|
||||
<DT><I>Cannot use pair hybrid with GPU neighbor builds</I>
|
||||
|
||||
<DD>Self-explanatory.
|
||||
<DD>See documentation for fix gpu.
|
||||
|
||||
<DT><I>Cannot use pair tail corrections with 2d simulations</I>
|
||||
|
||||
@ -1843,7 +1847,7 @@ does not exist.
|
||||
|
||||
<DD>Self-explanatory.
|
||||
|
||||
<DT><I>Could not find or initialize a specified accelerator device</I>
|
||||
<DT><I>Could not find/initialize a specified accelerator device</I>
|
||||
|
||||
<DD>Your GPU setup is invalid.
|
||||
|
||||
@ -2123,6 +2127,10 @@ model.
|
||||
used. Most likely, one or more atoms have been blown out of the
|
||||
simulation box to a great distance.
|
||||
|
||||
<DT><I>Double precision is not supported on this accelerator.</I>
|
||||
|
||||
<DD>In this case, you must compile the GPU library for single precision.
|
||||
|
||||
<DT><I>Dump cfg and fix not computed at compatible times</I>
|
||||
|
||||
<DD>The fix must produce per-atom quantities on timesteps that dump cfg
|
||||
@ -2355,6 +2363,10 @@ smaller simulation or on more processors.
|
||||
|
||||
<DD>Self-explanatory.
|
||||
|
||||
<DT><I>Fix gpu split must be positive for hybrid pair styles.</I>
|
||||
|
||||
<DD>See documentation for fix gpu.
|
||||
|
||||
<DT><I>Fix ID for compute atom/molecule does not exist</I>
|
||||
|
||||
<DD>Self-explanatory.
|
||||
@ -3227,6 +3239,11 @@ this fix.
|
||||
|
||||
<DD>This is the way the fix must be defined in your input script.
|
||||
|
||||
<DT><I>GPU library not compiled for this accelerator</I>
|
||||
|
||||
<DD>The GPU library was not built for your accelerator. Check the arch flag in
|
||||
lib/gpu.
|
||||
|
||||
<DT><I>Gmask function in equal-style variable formula</I>
|
||||
|
||||
<DD>Gmask is per-atom operation.
|
||||
@ -3509,7 +3526,7 @@ simulation box.
|
||||
|
||||
<DD>Eigensolve for rigid body was not sufficiently accurate.
|
||||
|
||||
<DT><I>Insufficient memory on accelerator (or no fix gpu)</I>
|
||||
<DT><I>Insufficient memory on accelerator. </I>
|
||||
|
||||
<DD>Self-explanatory.
|
||||
|
||||
@ -4587,10 +4604,6 @@ contain the same atom.
|
||||
<DD>Any rigid body defined by the fix rigid command must contain 2 or more
|
||||
atoms.
|
||||
|
||||
<DT><I>Out of memory on GPGPU</I>
|
||||
|
||||
<DD>You are attempting to run with too many atoms on the GPU.
|
||||
|
||||
<DT><I>Out of range atoms - cannot compute PPPM</I>
|
||||
|
||||
<DD>One or more atoms are attempting to map their charge to a PPPM grid
|
||||
|
||||
@ -505,6 +505,14 @@ the list.
|
||||
|
||||
|
||||
<DIV ALIGN=center><TABLE BORDER=1 >
|
||||
<TR><TD >pppm GPU single and double </TD><TD > Mike Brown (ORNL)</TD></TR>
|
||||
<TR><TD >pair_style lj/cut/expand </TD><TD > Inderaj Bains (NVIDIA)</TD></TR>
|
||||
<TR><TD >temperature accelerated dynamics (TAD) </TD><TD > Aidan Thompson (Sandia)</TD></TR>
|
||||
<TR><TD >pair reax/c and fix qeq/reax </TD><TD > Metin Aktulga (Purdue, now LBNL)</TD></TR>
|
||||
<TR><TD >DREIDING force field, pair_style hbond/dreiding, etc </TD><TD > Tod Pascal (CalTech)</TD></TR>
|
||||
<TR><TD >fix adapt and compute ti for thermodynamic integreation for free energies </TD><TD > Sai Jayaraman (Sandia)</TD></TR>
|
||||
<TR><TD >pair born and pair gauss </TD><TD > Sai Jayaraman (Sandia)</TD></TR>
|
||||
<TR><TD >stochastic rotation dynamics (SRD) via fix srd </TD><TD > Jemery Lechman (Sandia) and Pieter in 't Veld (BASF)</TD></TR>
|
||||
<TR><TD >ipp Perl script tool </TD><TD > Reese Jones (Sandia)</TD></TR>
|
||||
<TR><TD >eam_database and createatoms tools </TD><TD > Xiaowang Zhou (Sandia)</TD></TR>
|
||||
<TR><TD >electron force field (eFF) </TD><TD > Andres Jaramillo-Botero and Julius Su (Caltech)</TD></TR>
|
||||
|
||||
@ -490,7 +490,14 @@ the list.
|
||||
|
||||
:link(sjp,http://www.sandia.gov/~sjplimp)
|
||||
|
||||
pppm GPU single and double : Mike Brown (ORNL)
|
||||
pair_style lj/cut/expand : Inderaj Bains (NVIDIA)
|
||||
temperature accelerated dynamics (TAD) : Aidan Thompson (Sandia)
|
||||
pair reax/c and fix qeq/reax : Metin Aktulga (Purdue, now LBNL)
|
||||
DREIDING force field, pair_style hbond/dreiding, etc : Tod Pascal (CalTech)
|
||||
fix adapt and compute ti for thermodynamic integreation for free energies : Sai Jayaraman (Sandia)
|
||||
pair born and pair gauss : Sai Jayaraman (Sandia)
|
||||
stochastic rotation dynamics (SRD) via fix srd : Jemery Lechman (Sandia) and Pieter in 't Veld (BASF)
|
||||
ipp Perl script tool : Reese Jones (Sandia)
|
||||
eam_database and createatoms tools : Xiaowang Zhou (Sandia)
|
||||
electron force field (eFF) : Andres Jaramillo-Botero and Julius Su (Caltech)
|
||||
|
||||
@ -994,143 +994,130 @@ processing units (GPUs). We plan to add more over time. Currently,
|
||||
they only support NVIDIA GPU cards. To use them you need to install
|
||||
certain NVIDIA CUDA software on your system:
|
||||
</P>
|
||||
<UL><LI>Check if you have an NVIDIA card: cat /proc/driver/nvidia/cards/0
|
||||
<LI>Go to http://www.nvidia.com/object/cuda_get.html
|
||||
<LI>Install a driver and toolkit appropriate for your system (SDK is not necessary)
|
||||
<LI>Follow the instructions in README in lammps/lib/gpu to build the library.
|
||||
<LI>Run lammps/lib/gpu/nvc_get_devices to list supported devices and properties
|
||||
<UL><LI>Check if you have an NVIDIA card: cat /proc/driver/nvidia/cards/0 Go
|
||||
<LI>to http://www.nvidia.com/object/cuda_get.html Install a driver and
|
||||
<LI>toolkit appropriate for your system (SDK is not necessary) Follow the
|
||||
<LI>instructions in README in lammps/lib/gpu to build the library. Run
|
||||
<LI>lammps/lib/gpu/nvc_get_devices to list supported devices and
|
||||
<LI>properties
|
||||
</UL>
|
||||
<H4>GPU configuration
|
||||
</H4>
|
||||
<P>When using GPUs, you are restricted to one physical GPU per LAMMPS
|
||||
process. Multiple processes can share a single GPU and in many cases it
|
||||
will be more efficient to run with multiple processes per GPU. Any GPU
|
||||
accelerated style requires that <A HREF = "fix_gpu.html">fix gpu</A> be used in the
|
||||
input script to select and initialize the GPUs. The format for the fix
|
||||
is:
|
||||
process. Multiple processes can share a single GPU and in many cases
|
||||
it will be more efficient to run with multiple processes per GPU. Any
|
||||
GPU accelerated style requires that <A HREF = "fix_gpu.html">fix gpu</A> be used in
|
||||
the input script to select and initialize the GPUs. The format for the
|
||||
fix is:
|
||||
</P>
|
||||
<PRE>fix <I>name</I> all gpu <I>mode</I> <I>first</I> <I>last</I> <I>split</I>
|
||||
</PRE>
|
||||
<P>where <I>name</I> is the name for the fix. The gpu fix must be the first
|
||||
fix specified for a given run, otherwise the program will exit
|
||||
with an error. The gpu fix will not have any effect on runs
|
||||
that do not use GPU acceleration; there should be no problem
|
||||
with specifying the fix first in any input script.
|
||||
fix specified for a given run, otherwise the program will exit with an
|
||||
error. The gpu fix will not have any effect on runs that do not use
|
||||
GPU acceleration; there should be no problem with specifying the fix
|
||||
first in any input script.
|
||||
</P>
|
||||
<P><I>mode</I> can be either "force" or "force/neigh". In the former,
|
||||
neighbor list calculation is performed on the CPU using the
|
||||
standard LAMMPS routines. In the latter, the neighbor list
|
||||
calculation is performed on the GPU. The GPU neighbor list
|
||||
can be used for better performance, however, it
|
||||
should not be used with a triclinic box.
|
||||
<P><I>mode</I> can be either "force" or "force/neigh". In the former, neighbor
|
||||
list calculation is performed on the CPU using the standard LAMMPS
|
||||
routines. In the latter, the neighbor list calculation is performed on
|
||||
the GPU. The GPU neighbor list can be used for better performance,
|
||||
however, it cannot not be used with a triclinic box or with
|
||||
<A HREF = "pair_hybrid.html">hybrid</A> pair styles.
|
||||
</P>
|
||||
<P>There are cases when it might be more efficient to select the CPU for neighbor
|
||||
list builds. If a non-GPU enabled style requires a neighbor list, it will also
|
||||
be built using CPU routines. Redundant CPU and GPU neighbor list calculations
|
||||
will typically be less efficient. For <A HREF = "pair_hybrid.html">hybrid</A> pair
|
||||
styles, GPU calculated neighbor lists might be less efficient because
|
||||
no particles will be skipped in a given neighbor list.
|
||||
<P>There are cases when it might be more efficient to select the CPU for
|
||||
neighbor list builds. If a non-GPU enabled style requires a neighbor
|
||||
list, it will also be built using CPU routines. Redundant CPU and GPU
|
||||
neighbor list calculations will typically be less efficient.
|
||||
</P>
|
||||
<P><I>first</I> is the ID (as reported by lammps/lib/gpu/nvc_get_devices)
|
||||
of the first GPU that will be used on each node. <I>last</I> is the
|
||||
ID of the last GPU that will be used on each node. If you have
|
||||
only one GPU per node, <I>first</I> and <I>last</I> will typically both be
|
||||
0. Selecting a non-sequential set of GPU IDs (e.g. 0,1,3)
|
||||
is not currently supported.
|
||||
<P><I>first</I> is the ID (as reported by lammps/lib/gpu/nvc_get_devices) of
|
||||
the first GPU that will be used on each node. <I>last</I> is the ID of the
|
||||
last GPU that will be used on each node. If you have only one GPU per
|
||||
node, <I>first</I> and <I>last</I> will typically both be 0. Selecting a
|
||||
non-sequential set of GPU IDs (e.g. 0,1,3) is not currently supported.
|
||||
</P>
|
||||
<P><I>split</I> is the fraction of particles whose forces, torques,
|
||||
energies, and/or virials will be calculated on the GPU. This
|
||||
can be used to perform CPU and GPU force calculations
|
||||
simultaneously. If <I>split</I> is negative, the software will
|
||||
attempt to calculate the optimal fraction automatically
|
||||
every 25 timesteps based on CPU and GPU timings. Because the GPU speedups
|
||||
are dependent on the number of particles, automatic calculation of the
|
||||
split can be less efficient, but typically results in loop times
|
||||
within 20% of an optimal fixed split.
|
||||
<P><I>split</I> is the fraction of particles whose forces, torques, energies,
|
||||
and/or virials will be calculated on the GPU. This can be used to
|
||||
perform CPU and GPU force calculations simultaneously. If <I>split</I> is
|
||||
negative, the software will attempt to calculate the optimal fraction
|
||||
automatically every 25 timesteps based on CPU and GPU timings. Because
|
||||
the GPU speedups are dependent on the number of particles, automatic
|
||||
calculation of the split can be less efficient, but typically results
|
||||
in loop times within 20% of an optimal fixed split.
|
||||
</P>
|
||||
<P>If you have two GPUs per node, 8 CPU cores per node, and
|
||||
would like to run on 4 nodes with dynamic balancing of
|
||||
force calculation across CPU and GPU cores, the fix
|
||||
might be
|
||||
<P>If you have two GPUs per node, 8 CPU cores per node, and would like to
|
||||
run on 4 nodes with dynamic balancing of force calculation across CPU
|
||||
and GPU cores, the fix might be
|
||||
</P>
|
||||
<PRE>fix 0 all gpu force/neigh 0 1 -1
|
||||
</PRE>
|
||||
<P>with LAMMPS run on 32 processes. In this case, all
|
||||
CPU cores and GPU devices on the nodes would be utilized.
|
||||
Each GPU device would be shared by 4 CPU cores. The
|
||||
CPU cores would perform force calculations for some
|
||||
fraction of the particles at the same time the GPUs
|
||||
performed force calculation for the other particles.
|
||||
<P>with LAMMPS run on 32 processes. In this case, all CPU cores and GPU
|
||||
devices on the nodes would be utilized. Each GPU device would be
|
||||
shared by 4 CPU cores. The CPU cores would perform force calculations
|
||||
for some fraction of the particles at the same time the GPUs performed
|
||||
force calculation for the other particles.
|
||||
</P>
|
||||
<P>Because of the large number of cores on each GPU
|
||||
device, it might be more efficient to run on fewer
|
||||
processes per GPU when the number of particles per process
|
||||
is small (100's of particles); this can be necessary
|
||||
to keep the GPU cores busy.
|
||||
<P>Because of the large number of cores on each GPU device, it might be
|
||||
more efficient to run on fewer processes per GPU when the number of
|
||||
particles per process is small (100's of particles); this can be
|
||||
necessary to keep the GPU cores busy.
|
||||
</P>
|
||||
<H4>GPU input script
|
||||
</H4>
|
||||
<P>In order to use GPU acceleration in LAMMPS,
|
||||
<A HREF = "fix_gpu.html">fix_gpu</A>
|
||||
should be used in order to initialize and configure the
|
||||
GPUs for use. Additionally, GPU enabled styles must be
|
||||
selected in the input script. Currently,
|
||||
this is limited to a few <A HREF = "pair_style.html">pair styles</A>.
|
||||
Some GPU-enabled styles have additional restrictions
|
||||
listed in their documentation.
|
||||
<P>In order to use GPU acceleration in LAMMPS, <A HREF = "fix_gpu.html">fix_gpu</A>
|
||||
should be used in order to initialize and configure the GPUs for
|
||||
use. Additionally, GPU enabled styles must be selected in the input
|
||||
script. Currently, this is limited to a few <A HREF = "pair_style.html">pair
|
||||
styles</A> and PPPM. Some GPU-enabled styles have
|
||||
additional restrictions listed in their documentation.
|
||||
</P>
|
||||
<H4>GPU asynchronous pair computation
|
||||
</H4>
|
||||
<P>The GPU accelerated pair styles can be used to perform
|
||||
pair style force calculation on the GPU while other
|
||||
calculations are
|
||||
performed on the CPU. One method to do this is to specify
|
||||
a <I>split</I> in the gpu fix as described above. In this case,
|
||||
force calculation for the pair style will also be performed
|
||||
on the CPU.
|
||||
<P>The GPU accelerated pair styles can be used to perform pair style
|
||||
force calculation on the GPU while other calculations are performed on
|
||||
the CPU. One method to do this is to specify a <I>split</I> in the gpu fix
|
||||
as described above. In this case, force calculation for the pair
|
||||
style will also be performed on the CPU.
|
||||
</P>
|
||||
<P>When the CPU work in a GPU pair style has finished,
|
||||
the next force computation will begin, possibly before the
|
||||
GPU has finished. If <I>split</I> is 1.0 in the gpu fix, the next
|
||||
force computation will begin almost immediately. This can
|
||||
be used to run a <A HREF = "pair_hybrid.html">hybrid</A> GPU pair style at
|
||||
the same time as a hybrid CPU pair style. In this case, the
|
||||
GPU pair style should be first in the hybrid command in order to
|
||||
perform simultaneous calculations. This also
|
||||
allows <A HREF = "bond_style.html">bond</A>, <A HREF = "angle_style.html">angle</A>,
|
||||
<A HREF = "dihedral_style.html">dihedral</A>, <A HREF = "improper_style.html">improper</A>,
|
||||
and <A HREF = "kspace_style.html">long-range</A> force
|
||||
computations to be run simultaneously with the GPU pair style.
|
||||
Once all CPU force computations have completed, the gpu fix
|
||||
will block until the GPU has finished all work before continuing
|
||||
the run.
|
||||
<P>When the CPU work in a GPU pair style has finished, the next force
|
||||
computation will begin, possibly before the GPU has finished. If
|
||||
<I>split</I> is 1.0 in the gpu fix, the next force computation will begin
|
||||
almost immediately. This can be used to run a
|
||||
<A HREF = "pair_hybrid.html">hybrid</A> GPU pair style at the same time as a hybrid
|
||||
CPU pair style. In this case, the GPU pair style should be first in
|
||||
the hybrid command in order to perform simultaneous calculations. This
|
||||
also allows <A HREF = "bond_style.html">bond</A>, <A HREF = "angle_style.html">angle</A>,
|
||||
<A HREF = "dihedral_style.html">dihedral</A>, <A HREF = "improper_style.html">improper</A>, and
|
||||
<A HREF = "kspace_style.html">long-range</A> force computations to be run
|
||||
simultaneously with the GPU pair style. Once all CPU force
|
||||
computations have completed, the gpu fix will block until the GPU has
|
||||
finished all work before continuing the run.
|
||||
</P>
|
||||
<H4>GPU timing
|
||||
</H4>
|
||||
<P>GPU accelerated pair styles can perform computations asynchronously
|
||||
with CPU computations. The "Pair" time reported by LAMMPS
|
||||
will be the maximum of the time required to complete the CPU
|
||||
pair style computations and the time required to complete the GPU
|
||||
pair style computations. Any time spent for GPU-enabled pair styles
|
||||
for computations that run simultaneously with <A HREF = "bond_style.html">bond</A>,
|
||||
<A HREF = "angle_style.html">angle</A>, <A HREF = "dihedral_style.html">dihedral</A>,
|
||||
<A HREF = "improper_style.html">improper</A>, and <A HREF = "kspace_style.html">long-range</A> calculations
|
||||
will not be included in the "Pair" time.
|
||||
with CPU computations. The "Pair" time reported by LAMMPS will be the
|
||||
maximum of the time required to complete the CPU pair style
|
||||
computations and the time required to complete the GPU pair style
|
||||
computations. Any time spent for GPU-enabled pair styles for
|
||||
computations that run simultaneously with <A HREF = "bond_style.html">bond</A>,
|
||||
<A HREF = "angle_style.html">angle</A>, <A HREF = "dihedral_style.html">dihedral</A>,
|
||||
<A HREF = "improper_style.html">improper</A>, and <A HREF = "kspace_style.html">long-range</A>
|
||||
calculations will not be included in the "Pair" time.
|
||||
</P>
|
||||
<P>When <I>mode</I> for the gpu fix is force/neigh,
|
||||
the time for neighbor list calculations on the GPU will be added
|
||||
into the "Pair" time, not the "Neigh" time. A breakdown of the
|
||||
times required for various tasks on the GPU (data copy, neighbor
|
||||
calculations, force computations, etc.) are output only
|
||||
with the LAMMPS screen output at the end of each run. These timings represent
|
||||
total time spent on the GPU for each routine, regardless of asynchronous
|
||||
CPU calculations.
|
||||
<P>When <I>mode</I> for the gpu fix is force/neigh, the time for neighbor list
|
||||
calculations on the GPU will be added into the "Pair" time, not the
|
||||
"Neigh" time. A breakdown of the times required for various tasks on
|
||||
the GPU (data copy, neighbor calculations, force computations, etc.)
|
||||
are output only with the LAMMPS screen output at the end of each
|
||||
run. These timings represent total time spent on the GPU for each
|
||||
routine, regardless of asynchronous CPU calculations.
|
||||
</P>
|
||||
<H4>GPU single vs double precision
|
||||
</H4>
|
||||
<P>See the lammps/lib/gpu/README file for instructions on how to build
|
||||
the LAMMPS gpu library for single, mixed, and double precision. The latter
|
||||
requires that your GPU card supports double precision.
|
||||
<P>See the lammps/lib/gpu/README file for instructions on how to build
|
||||
the LAMMPS gpu library for single, mixed, and double precision. The
|
||||
latter requires that your GPU card supports double precision.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
|
||||
@ -984,141 +984,130 @@ processing units (GPUs). We plan to add more over time. Currently,
|
||||
they only support NVIDIA GPU cards. To use them you need to install
|
||||
certain NVIDIA CUDA software on your system:
|
||||
|
||||
Check if you have an NVIDIA card: cat /proc/driver/nvidia/cards/0
|
||||
Go to http://www.nvidia.com/object/cuda_get.html
|
||||
Install a driver and toolkit appropriate for your system (SDK is not necessary)
|
||||
Follow the instructions in README in lammps/lib/gpu to build the library.
|
||||
Run lammps/lib/gpu/nvc_get_devices to list supported devices and properties :ul
|
||||
Check if you have an NVIDIA card: cat /proc/driver/nvidia/cards/0 Go
|
||||
to http://www.nvidia.com/object/cuda_get.html Install a driver and
|
||||
toolkit appropriate for your system (SDK is not necessary) Follow the
|
||||
instructions in README in lammps/lib/gpu to build the library. Run
|
||||
lammps/lib/gpu/nvc_get_devices to list supported devices and
|
||||
properties :ul
|
||||
|
||||
GPU configuration :h4
|
||||
|
||||
When using GPUs, you are restricted to one physical GPU per LAMMPS
|
||||
process. Multiple processes can share a single GPU and in many cases it
|
||||
will be more efficient to run with multiple processes per GPU. Any GPU
|
||||
accelerated style requires that "fix gpu"_fix_gpu.html be used in the
|
||||
input script to select and initialize the GPUs. The format for the fix
|
||||
is:
|
||||
process. Multiple processes can share a single GPU and in many cases
|
||||
it will be more efficient to run with multiple processes per GPU. Any
|
||||
GPU accelerated style requires that "fix gpu"_fix_gpu.html be used in
|
||||
the input script to select and initialize the GPUs. The format for the
|
||||
fix is:
|
||||
|
||||
fix {name} all gpu {mode} {first} {last} {split} :pre
|
||||
|
||||
where {name} is the name for the fix. The gpu fix must be the first
|
||||
fix specified for a given run, otherwise the program will exit
|
||||
with an error. The gpu fix will not have any effect on runs
|
||||
that do not use GPU acceleration; there should be no problem
|
||||
with specifying the fix first in any input script.
|
||||
fix specified for a given run, otherwise the program will exit with an
|
||||
error. The gpu fix will not have any effect on runs that do not use
|
||||
GPU acceleration; there should be no problem with specifying the fix
|
||||
first in any input script.
|
||||
|
||||
{mode} can be either "force" or "force/neigh". In the former,
|
||||
neighbor list calculation is performed on the CPU using the
|
||||
standard LAMMPS routines. In the latter, the neighbor list
|
||||
calculation is performed on the GPU. The GPU neighbor list
|
||||
can be used for better performance, however, it
|
||||
cannot not be used with a triclinic box or with "hybrid"_pair_hybrid.html
|
||||
pair styles.
|
||||
{mode} can be either "force" or "force/neigh". In the former, neighbor
|
||||
list calculation is performed on the CPU using the standard LAMMPS
|
||||
routines. In the latter, the neighbor list calculation is performed on
|
||||
the GPU. The GPU neighbor list can be used for better performance,
|
||||
however, it cannot not be used with a triclinic box or with
|
||||
"hybrid"_pair_hybrid.html pair styles.
|
||||
|
||||
There are cases when it might be more efficient to select the CPU for neighbor
|
||||
list builds. If a non-GPU enabled style requires a neighbor list, it will also
|
||||
be built using CPU routines. Redundant CPU and GPU neighbor list calculations
|
||||
will typically be less efficient.
|
||||
There are cases when it might be more efficient to select the CPU for
|
||||
neighbor list builds. If a non-GPU enabled style requires a neighbor
|
||||
list, it will also be built using CPU routines. Redundant CPU and GPU
|
||||
neighbor list calculations will typically be less efficient.
|
||||
|
||||
{first} is the ID (as reported by lammps/lib/gpu/nvc_get_devices)
|
||||
of the first GPU that will be used on each node. {last} is the
|
||||
ID of the last GPU that will be used on each node. If you have
|
||||
only one GPU per node, {first} and {last} will typically both be
|
||||
0. Selecting a non-sequential set of GPU IDs (e.g. 0,1,3)
|
||||
is not currently supported.
|
||||
{first} is the ID (as reported by lammps/lib/gpu/nvc_get_devices) of
|
||||
the first GPU that will be used on each node. {last} is the ID of the
|
||||
last GPU that will be used on each node. If you have only one GPU per
|
||||
node, {first} and {last} will typically both be 0. Selecting a
|
||||
non-sequential set of GPU IDs (e.g. 0,1,3) is not currently supported.
|
||||
|
||||
{split} is the fraction of particles whose forces, torques,
|
||||
energies, and/or virials will be calculated on the GPU. This
|
||||
can be used to perform CPU and GPU force calculations
|
||||
simultaneously. If {split} is negative, the software will
|
||||
attempt to calculate the optimal fraction automatically
|
||||
every 25 timesteps based on CPU and GPU timings. Because the GPU speedups
|
||||
are dependent on the number of particles, automatic calculation of the
|
||||
split can be less efficient, but typically results in loop times
|
||||
within 20% of an optimal fixed split.
|
||||
{split} is the fraction of particles whose forces, torques, energies,
|
||||
and/or virials will be calculated on the GPU. This can be used to
|
||||
perform CPU and GPU force calculations simultaneously. If {split} is
|
||||
negative, the software will attempt to calculate the optimal fraction
|
||||
automatically every 25 timesteps based on CPU and GPU timings. Because
|
||||
the GPU speedups are dependent on the number of particles, automatic
|
||||
calculation of the split can be less efficient, but typically results
|
||||
in loop times within 20% of an optimal fixed split.
|
||||
|
||||
If you have two GPUs per node, 8 CPU cores per node, and
|
||||
would like to run on 4 nodes with dynamic balancing of
|
||||
force calculation across CPU and GPU cores, the fix
|
||||
might be
|
||||
If you have two GPUs per node, 8 CPU cores per node, and would like to
|
||||
run on 4 nodes with dynamic balancing of force calculation across CPU
|
||||
and GPU cores, the fix might be
|
||||
|
||||
fix 0 all gpu force/neigh 0 1 -1 :pre
|
||||
|
||||
with LAMMPS run on 32 processes. In this case, all
|
||||
CPU cores and GPU devices on the nodes would be utilized.
|
||||
Each GPU device would be shared by 4 CPU cores. The
|
||||
CPU cores would perform force calculations for some
|
||||
fraction of the particles at the same time the GPUs
|
||||
performed force calculation for the other particles.
|
||||
with LAMMPS run on 32 processes. In this case, all CPU cores and GPU
|
||||
devices on the nodes would be utilized. Each GPU device would be
|
||||
shared by 4 CPU cores. The CPU cores would perform force calculations
|
||||
for some fraction of the particles at the same time the GPUs performed
|
||||
force calculation for the other particles.
|
||||
|
||||
Because of the large number of cores on each GPU
|
||||
device, it might be more efficient to run on fewer
|
||||
processes per GPU when the number of particles per process
|
||||
is small (100's of particles); this can be necessary
|
||||
to keep the GPU cores busy.
|
||||
Because of the large number of cores on each GPU device, it might be
|
||||
more efficient to run on fewer processes per GPU when the number of
|
||||
particles per process is small (100's of particles); this can be
|
||||
necessary to keep the GPU cores busy.
|
||||
|
||||
GPU input script :h4
|
||||
|
||||
In order to use GPU acceleration in LAMMPS,
|
||||
"fix_gpu"_fix_gpu.html
|
||||
should be used in order to initialize and configure the
|
||||
GPUs for use. Additionally, GPU enabled styles must be
|
||||
selected in the input script. Currently, this is limited
|
||||
to a few "pair styles"_pair_style.html and PPPM.
|
||||
Some GPU-enabled styles have additional restrictions
|
||||
listed in their documentation.
|
||||
In order to use GPU acceleration in LAMMPS, "fix_gpu"_fix_gpu.html
|
||||
should be used in order to initialize and configure the GPUs for
|
||||
use. Additionally, GPU enabled styles must be selected in the input
|
||||
script. Currently, this is limited to a few "pair
|
||||
styles"_pair_style.html and PPPM. Some GPU-enabled styles have
|
||||
additional restrictions listed in their documentation.
|
||||
|
||||
GPU asynchronous pair computation :h4
|
||||
|
||||
The GPU accelerated pair styles can be used to perform
|
||||
pair style force calculation on the GPU while other
|
||||
calculations are performed on the CPU. One method to do this
|
||||
is to specify a {split} in the gpu fix as described above.
|
||||
In this case, force calculation for the pair style will also
|
||||
be performed on the CPU.
|
||||
The GPU accelerated pair styles can be used to perform pair style
|
||||
force calculation on the GPU while other calculations are performed on
|
||||
the CPU. One method to do this is to specify a {split} in the gpu fix
|
||||
as described above. In this case, force calculation for the pair
|
||||
style will also be performed on the CPU.
|
||||
|
||||
When the CPU work in a GPU pair style has finished,
|
||||
the next force computation will begin, possibly before the
|
||||
GPU has finished. If {split} is 1.0 in the gpu fix, the next
|
||||
force computation will begin almost immediately. This can
|
||||
be used to run a "hybrid"_pair_hybrid.html GPU pair style at
|
||||
the same time as a hybrid CPU pair style. In this case, the
|
||||
GPU pair style should be first in the hybrid command in order to
|
||||
perform simultaneous calculations. This also
|
||||
allows "bond"_bond_style.html, "angle"_angle_style.html,
|
||||
"dihedral"_dihedral_style.html, "improper"_improper_style.html,
|
||||
and "long-range"_kspace_style.html force
|
||||
computations to be run simultaneously with the GPU pair style.
|
||||
Once all CPU force computations have completed, the gpu fix
|
||||
will block until the GPU has finished all work before continuing
|
||||
the run.
|
||||
When the CPU work in a GPU pair style has finished, the next force
|
||||
computation will begin, possibly before the GPU has finished. If
|
||||
{split} is 1.0 in the gpu fix, the next force computation will begin
|
||||
almost immediately. This can be used to run a
|
||||
"hybrid"_pair_hybrid.html GPU pair style at the same time as a hybrid
|
||||
CPU pair style. In this case, the GPU pair style should be first in
|
||||
the hybrid command in order to perform simultaneous calculations. This
|
||||
also allows "bond"_bond_style.html, "angle"_angle_style.html,
|
||||
"dihedral"_dihedral_style.html, "improper"_improper_style.html, and
|
||||
"long-range"_kspace_style.html force computations to be run
|
||||
simultaneously with the GPU pair style. Once all CPU force
|
||||
computations have completed, the gpu fix will block until the GPU has
|
||||
finished all work before continuing the run.
|
||||
|
||||
GPU timing :h4
|
||||
|
||||
GPU accelerated pair styles can perform computations asynchronously
|
||||
with CPU computations. The "Pair" time reported by LAMMPS
|
||||
will be the maximum of the time required to complete the CPU
|
||||
pair style computations and the time required to complete the GPU
|
||||
pair style computations. Any time spent for GPU-enabled pair styles
|
||||
for computations that run simultaneously with "bond"_bond_style.html,
|
||||
"angle"_angle_style.html, "dihedral"_dihedral_style.html,
|
||||
"improper"_improper_style.html, and "long-range"_kspace_style.html calculations
|
||||
will not be included in the "Pair" time.
|
||||
with CPU computations. The "Pair" time reported by LAMMPS will be the
|
||||
maximum of the time required to complete the CPU pair style
|
||||
computations and the time required to complete the GPU pair style
|
||||
computations. Any time spent for GPU-enabled pair styles for
|
||||
computations that run simultaneously with "bond"_bond_style.html,
|
||||
"angle"_angle_style.html, "dihedral"_dihedral_style.html,
|
||||
"improper"_improper_style.html, and "long-range"_kspace_style.html
|
||||
calculations will not be included in the "Pair" time.
|
||||
|
||||
When {mode} for the gpu fix is force/neigh,
|
||||
the time for neighbor list calculations on the GPU will be added
|
||||
into the "Pair" time, not the "Neigh" time. A breakdown of the
|
||||
times required for various tasks on the GPU (data copy, neighbor
|
||||
calculations, force computations, etc.) are output only
|
||||
with the LAMMPS screen output at the end of each run. These timings represent
|
||||
total time spent on the GPU for each routine, regardless of asynchronous
|
||||
CPU calculations.
|
||||
When {mode} for the gpu fix is force/neigh, the time for neighbor list
|
||||
calculations on the GPU will be added into the "Pair" time, not the
|
||||
"Neigh" time. A breakdown of the times required for various tasks on
|
||||
the GPU (data copy, neighbor calculations, force computations, etc.)
|
||||
are output only with the LAMMPS screen output at the end of each
|
||||
run. These timings represent total time spent on the GPU for each
|
||||
routine, regardless of asynchronous CPU calculations.
|
||||
|
||||
GPU single vs double precision :h4
|
||||
|
||||
See the lammps/lib/gpu/README file for instructions on how to build
|
||||
the LAMMPS gpu library for single, mixed, and double precision. The latter
|
||||
requires that your GPU card supports double precision.
|
||||
See the lammps/lib/gpu/README file for instructions on how to build
|
||||
the LAMMPS gpu library for single, mixed, and double precision. The
|
||||
latter requires that your GPU card supports double precision.
|
||||
|
||||
:line
|
||||
|
||||
|
||||
@ -13,16 +13,29 @@
|
||||
</H3>
|
||||
<P><B>Syntax:</B>
|
||||
</P>
|
||||
<PRE>compute ID group-ID temp/asphere bias-ID
|
||||
<PRE>compute ID group-ID temp/asphere keyword value ...
|
||||
</PRE>
|
||||
<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
|
||||
<LI>temp/asphere = style name of this compute command
|
||||
<LI>bias-ID = ID of a temperature compute that removes a velocity bias (optional)
|
||||
<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
|
||||
|
||||
<LI>temp/asphere = style name of this compute command
|
||||
|
||||
<LI>zero or more keyword/value pairs may be appended
|
||||
|
||||
<LI>keyword = <I>bias</I> or <I>dof</I>
|
||||
|
||||
<PRE> <I>bias</I> value = bias-ID<I>uniform</I> or <I>gaussian</I>
|
||||
bias-ID = ID of a temperature compute that removes a velocity bias
|
||||
<I>dof</I> value = <I>all</I> or <I>rotate</I>
|
||||
all = compute temperature of translational and rotational degrees of freedom
|
||||
rotate = compute temperature of just rotational degrees of freedom
|
||||
</PRE>
|
||||
|
||||
</UL>
|
||||
<P><B>Examples:</B>
|
||||
</P>
|
||||
<PRE>compute 1 all temp/asphere
|
||||
compute myTemp mobile temp/asphere tempCOM
|
||||
compute myTemp mobile temp/asphere bias tempCOM
|
||||
compute myTemp mobile temp/asphere dof rotate
|
||||
</PRE>
|
||||
<P><B>Description:</B>
|
||||
</P>
|
||||
@ -75,15 +88,6 @@ vector are ordered xx, yy, zz, xy, xz, yz.
|
||||
constant for the duration of the run; use the <I>dynamic</I> option of the
|
||||
<A HREF = "compute_modify.html">compute_modify</A> command if this is not the case.
|
||||
</P>
|
||||
<P>If a <I>bias-ID</I> is specified it must be the ID of a temperature compute
|
||||
that removes a "bias" velocity from each atom. This allows compute
|
||||
temp/sphere to compute its thermal temperature after the translational
|
||||
kinetic energy components have been altered in a prescribed way,
|
||||
e.g. to remove a velocity profile. Thermostats that use this compute
|
||||
will work with this bias term. See the doc pages for individual
|
||||
computes that calculate a temperature and the doc pages for fixes that
|
||||
perform thermostatting for more details.
|
||||
</P>
|
||||
<P>This compute subtracts out translational degrees-of-freedom due to
|
||||
fixes that constrain molecular motion, such as <A HREF = "fix_shake.html">fix
|
||||
shake</A> and <A HREF = "fix_rigid.html">fix rigid</A>. This means the
|
||||
@ -96,6 +100,26 @@ be altered using the <I>extra</I> option of the
|
||||
discussion of different ways to compute temperature and perform
|
||||
thermostatting.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The keyword/value option pairs are used in the following ways.
|
||||
</P>
|
||||
<P>For the <I>bias</I> keyword, <I>bias-ID</I> refers to the ID of a temperature
|
||||
compute that removes a "bias" velocity from each atom. This allows
|
||||
compute temp/sphere to compute its thermal temperature after the
|
||||
translational kinetic energy components have been altered in a
|
||||
prescribed way, e.g. to remove a velocity profile. Thermostats that
|
||||
use this compute will work with this bias term. See the doc pages for
|
||||
individual computes that calculate a temperature and the doc pages for
|
||||
fixes that perform thermostatting for more details.
|
||||
</P>
|
||||
<P>For the <I>dof</I> keyword, a setting of <I>all</I> calculates a temperature
|
||||
that includes both translational and rotational degrees of freedom. A
|
||||
setting of <I>rotate</I> calculates a temperature that includes only
|
||||
rotational degrees of freedom.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P><B>Output info:</B>
|
||||
</P>
|
||||
<P>This compute calculates a global scalar (the temperature) and a global
|
||||
|
||||
@ -10,16 +10,24 @@ compute temp/asphere command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
compute ID group-ID temp/asphere bias-ID :pre
|
||||
compute ID group-ID temp/asphere keyword value ... :pre
|
||||
|
||||
ID, group-ID are documented in "compute"_compute.html command
|
||||
temp/asphere = style name of this compute command
|
||||
bias-ID = ID of a temperature compute that removes a velocity bias (optional) :ul
|
||||
ID, group-ID are documented in "compute"_compute.html command :ulb,l
|
||||
temp/asphere = style name of this compute command :l
|
||||
zero or more keyword/value pairs may be appended :l
|
||||
keyword = {bias} or {dof} :l
|
||||
{bias} value = bias-ID{uniform} or {gaussian}
|
||||
bias-ID = ID of a temperature compute that removes a velocity bias
|
||||
{dof} value = {all} or {rotate}
|
||||
all = compute temperature of translational and rotational degrees of freedom
|
||||
rotate = compute temperature of just rotational degrees of freedom :pre
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
compute 1 all temp/asphere
|
||||
compute myTemp mobile temp/asphere tempCOM :pre
|
||||
compute myTemp mobile temp/asphere bias tempCOM
|
||||
compute myTemp mobile temp/asphere dof rotate :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
@ -72,15 +80,6 @@ The number of atoms contributing to the temperature is assumed to be
|
||||
constant for the duration of the run; use the {dynamic} option of the
|
||||
"compute_modify"_compute_modify.html command if this is not the case.
|
||||
|
||||
If a {bias-ID} is specified it must be the ID of a temperature compute
|
||||
that removes a "bias" velocity from each atom. This allows compute
|
||||
temp/sphere to compute its thermal temperature after the translational
|
||||
kinetic energy components have been altered in a prescribed way,
|
||||
e.g. to remove a velocity profile. Thermostats that use this compute
|
||||
will work with this bias term. See the doc pages for individual
|
||||
computes that calculate a temperature and the doc pages for fixes that
|
||||
perform thermostatting for more details.
|
||||
|
||||
This compute subtracts out translational degrees-of-freedom due to
|
||||
fixes that constrain molecular motion, such as "fix
|
||||
shake"_fix_shake.html and "fix rigid"_fix_rigid.html. This means the
|
||||
@ -93,6 +92,26 @@ See "this howto section"_Section_howto.html#4_16 of the manual for a
|
||||
discussion of different ways to compute temperature and perform
|
||||
thermostatting.
|
||||
|
||||
:line
|
||||
|
||||
The keyword/value option pairs are used in the following ways.
|
||||
|
||||
For the {bias} keyword, {bias-ID} refers to the ID of a temperature
|
||||
compute that removes a "bias" velocity from each atom. This allows
|
||||
compute temp/sphere to compute its thermal temperature after the
|
||||
translational kinetic energy components have been altered in a
|
||||
prescribed way, e.g. to remove a velocity profile. Thermostats that
|
||||
use this compute will work with this bias term. See the doc pages for
|
||||
individual computes that calculate a temperature and the doc pages for
|
||||
fixes that perform thermostatting for more details.
|
||||
|
||||
For the {dof} keyword, a setting of {all} calculates a temperature
|
||||
that includes both translational and rotational degrees of freedom. A
|
||||
setting of {rotate} calculates a temperature that includes only
|
||||
rotational degrees of freedom.
|
||||
|
||||
:line
|
||||
|
||||
[Output info:]
|
||||
|
||||
This compute calculates a global scalar (the temperature) and a global
|
||||
|
||||
@ -13,16 +13,29 @@
|
||||
</H3>
|
||||
<P><B>Syntax:</B>
|
||||
</P>
|
||||
<PRE>compute ID group-ID temp/sphere bias-ID
|
||||
<PRE>compute ID group-ID temp/sphere keyword value ...
|
||||
</PRE>
|
||||
<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
|
||||
<LI>temp/sphere = style name of this compute command
|
||||
<LI>bias-ID = ID of a temperature compute that removes a velocity bias (optional)
|
||||
<UL><LI>ID, group-ID are documented in <A HREF = "compute.html">compute</A> command
|
||||
|
||||
<LI>temp/sphere = style name of this compute command
|
||||
|
||||
<LI>zero or more keyword/value pairs may be appended
|
||||
|
||||
<LI>keyword = <I>bias</I> or <I>dof</I>
|
||||
|
||||
<PRE> <I>bias</I> value = bias-ID<I>uniform</I> or <I>gaussian</I>
|
||||
bias-ID = ID of a temperature compute that removes a velocity bias
|
||||
<I>dof</I> value = <I>all</I> or <I>rotate</I>
|
||||
all = compute temperature of translational and rotational degrees of freedom
|
||||
rotate = compute temperature of just rotational degrees of freedom
|
||||
</PRE>
|
||||
|
||||
</UL>
|
||||
<P><B>Examples:</B>
|
||||
</P>
|
||||
<PRE>compute 1 all temp/sphere
|
||||
compute myTemp mobile temp/sphere tempCOM
|
||||
compute myTemp mobile temp/sphere bias tempCOM
|
||||
compute myTemp mobile temp/sphere dof rotate
|
||||
</PRE>
|
||||
<P><B>Description:</B>
|
||||
</P>
|
||||
@ -66,15 +79,6 @@ the vector are ordered xx, yy, zz, xy, xz, yz.
|
||||
constant for the duration of the run; use the <I>dynamic</I> option of the
|
||||
<A HREF = "compute_modify.html">compute_modify</A> command if this is not the case.
|
||||
</P>
|
||||
<P>If a <I>bias-ID</I> is specified it must be the ID of a temperature compute
|
||||
that removes a "bias" velocity from each atom. This allows compute
|
||||
temp/sphere to compute its thermal temperature after the translational
|
||||
kinetic energy components have been altered in a prescribed way,
|
||||
e.g. to remove a velocity profile. Thermostats that use this compute
|
||||
will work with this bias term. See the doc pages for individual
|
||||
computes that calculate a temperature and the doc pages for fixes that
|
||||
perform thermostatting for more details.
|
||||
</P>
|
||||
<P>This compute subtracts out translational degrees-of-freedom due to
|
||||
fixes that constrain molecular motion, such as <A HREF = "fix_shake.html">fix
|
||||
shake</A> and <A HREF = "fix_rigid.html">fix rigid</A>. This means the
|
||||
@ -87,6 +91,26 @@ be altered using the <I>extra</I> option of the
|
||||
discussion of different ways to compute temperature and perform
|
||||
thermostatting.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The keyword/value option pairs are used in the following ways.
|
||||
</P>
|
||||
<P>For the <I>bias</I> keyword, <I>bias-ID</I> refers to the ID of a temperature
|
||||
compute that removes a "bias" velocity from each atom. This allows
|
||||
compute temp/sphere to compute its thermal temperature after the
|
||||
translational kinetic energy components have been altered in a
|
||||
prescribed way, e.g. to remove a velocity profile. Thermostats that
|
||||
use this compute will work with this bias term. See the doc pages for
|
||||
individual computes that calculate a temperature and the doc pages for
|
||||
fixes that perform thermostatting for more details.
|
||||
</P>
|
||||
<P>For the <I>dof</I> keyword, a setting of <I>all</I> calculates a temperature
|
||||
that includes both translational and rotational degrees of freedom. A
|
||||
setting of <I>rotate</I> calculates a temperature that includes only
|
||||
rotational degrees of freedom.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P><B>Output info:</B>
|
||||
</P>
|
||||
<P>This compute calculates a global scalar (the temperature) and a global
|
||||
@ -116,6 +140,8 @@ particles with radius = 0.0.
|
||||
<P><A HREF = "compute_temp.html">compute temp</A>, <A HREF = "compute_temp.html">compute
|
||||
temp/asphere</A>
|
||||
</P>
|
||||
<P><B>Default:</B> none
|
||||
<P><B>Default:</B>
|
||||
</P>
|
||||
<P>The option defaults are no bias and dof = all.
|
||||
</P>
|
||||
</HTML>
|
||||
|
||||
@ -10,16 +10,24 @@ compute temp/sphere command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
compute ID group-ID temp/sphere bias-ID :pre
|
||||
compute ID group-ID temp/sphere keyword value ... :pre
|
||||
|
||||
ID, group-ID are documented in "compute"_compute.html command
|
||||
temp/sphere = style name of this compute command
|
||||
bias-ID = ID of a temperature compute that removes a velocity bias (optional) :ul
|
||||
ID, group-ID are documented in "compute"_compute.html command :ulb,l
|
||||
temp/sphere = style name of this compute command :l
|
||||
zero or more keyword/value pairs may be appended :l
|
||||
keyword = {bias} or {dof} :l
|
||||
{bias} value = bias-ID{uniform} or {gaussian}
|
||||
bias-ID = ID of a temperature compute that removes a velocity bias
|
||||
{dof} value = {all} or {rotate}
|
||||
all = compute temperature of translational and rotational degrees of freedom
|
||||
rotate = compute temperature of just rotational degrees of freedom :pre
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
compute 1 all temp/sphere
|
||||
compute myTemp mobile temp/sphere tempCOM :pre
|
||||
compute myTemp mobile temp/sphere bias tempCOM
|
||||
compute myTemp mobile temp/sphere dof rotate :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
@ -63,15 +71,6 @@ The number of atoms contributing to the temperature is assumed to be
|
||||
constant for the duration of the run; use the {dynamic} option of the
|
||||
"compute_modify"_compute_modify.html command if this is not the case.
|
||||
|
||||
If a {bias-ID} is specified it must be the ID of a temperature compute
|
||||
that removes a "bias" velocity from each atom. This allows compute
|
||||
temp/sphere to compute its thermal temperature after the translational
|
||||
kinetic energy components have been altered in a prescribed way,
|
||||
e.g. to remove a velocity profile. Thermostats that use this compute
|
||||
will work with this bias term. See the doc pages for individual
|
||||
computes that calculate a temperature and the doc pages for fixes that
|
||||
perform thermostatting for more details.
|
||||
|
||||
This compute subtracts out translational degrees-of-freedom due to
|
||||
fixes that constrain molecular motion, such as "fix
|
||||
shake"_fix_shake.html and "fix rigid"_fix_rigid.html. This means the
|
||||
@ -84,6 +83,26 @@ See "this howto section"_Section_howto.html#4_16 of the manual for a
|
||||
discussion of different ways to compute temperature and perform
|
||||
thermostatting.
|
||||
|
||||
:line
|
||||
|
||||
The keyword/value option pairs are used in the following ways.
|
||||
|
||||
For the {bias} keyword, {bias-ID} refers to the ID of a temperature
|
||||
compute that removes a "bias" velocity from each atom. This allows
|
||||
compute temp/sphere to compute its thermal temperature after the
|
||||
translational kinetic energy components have been altered in a
|
||||
prescribed way, e.g. to remove a velocity profile. Thermostats that
|
||||
use this compute will work with this bias term. See the doc pages for
|
||||
individual computes that calculate a temperature and the doc pages for
|
||||
fixes that perform thermostatting for more details.
|
||||
|
||||
For the {dof} keyword, a setting of {all} calculates a temperature
|
||||
that includes both translational and rotational degrees of freedom. A
|
||||
setting of {rotate} calculates a temperature that includes only
|
||||
rotational degrees of freedom.
|
||||
|
||||
:line
|
||||
|
||||
[Output info:]
|
||||
|
||||
This compute calculates a global scalar (the temperature) and a global
|
||||
@ -113,4 +132,6 @@ particles with radius = 0.0.
|
||||
"compute temp"_compute_temp.html, "compute
|
||||
temp/asphere"_compute_temp.html
|
||||
|
||||
[Default:] none
|
||||
[Default:]
|
||||
|
||||
The option defaults are no bias and dof = all.
|
||||
|
||||
@ -47,7 +47,8 @@
|
||||
</PRE>
|
||||
<PRE> <I>custom</I> args = list of atom attributes
|
||||
possible attributes = id, mol, type, mass,
|
||||
x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz,
|
||||
x, y, z, xs, ys, zs, xu, yu, zu,
|
||||
xsu, ysu, zsu, ix, iy, iz,
|
||||
vx, vy, vz, fx, fy, fz,
|
||||
q, mux, muy, muz, mu,
|
||||
radius, omegax, omegay, omegaz,
|
||||
@ -62,6 +63,7 @@
|
||||
x,y,z = unscaled atom coordinates
|
||||
xs,ys,zs = scaled atom coordinates
|
||||
xu,yu,zu = unwrapped atom coordinates
|
||||
xsu,ysu,zsu = scaled unwrapped atom coordinates
|
||||
ix,iy,iz = box image that the atom is in
|
||||
vx,vy,vz = atom velocities
|
||||
fx,fy,fz = forces on atoms
|
||||
@ -228,14 +230,23 @@ extended CFG format files, as used by the
|
||||
package. Since the extended CFG format uses a single snapshot of the
|
||||
system per file, a wildcard "*" must be included in the filename, as
|
||||
discussed below. The list of atom attributes for style <I>cfg</I> must
|
||||
begin with "id type xs ys zs", since these quantities are needed to
|
||||
begin with either "id type xs ys zs" or "id type xsu ysu zsu" or
|
||||
since these quantities are needed to
|
||||
write the CFG files in the appropriate format (though the "id" and
|
||||
"type" fields do not appear explicitly in the file). Any remaining
|
||||
attributes will be stored as "auxiliary properties" in the CFG files.
|
||||
Note that you will typically want to use the <A HREF = "dump_modify.html">dump_modify
|
||||
element</A> command with CFG-formatted files, to
|
||||
associate element names with atom types, so that AtomEye can render
|
||||
atoms appropriately.
|
||||
atoms appropriately. When unwrapped coordinates <I>xsu</I>, <I>ysu</I>, and <I>zsu</I>
|
||||
are requested, the nominal AtomEye periodic cell dimensions are expanded
|
||||
by a large factor UNWRAPEXPAND = 10.0, which ensures atoms that are
|
||||
displayed correctly for up to UNWRAPEXPAND/2 periodic boundary crossings
|
||||
in any direction.
|
||||
Beyond this, AtomEye will rewrap the unwrapped coordinates.
|
||||
The expansion causes the atoms to be drawn farther
|
||||
away from the viewer, but it is easy to zoom the atoms closer, and
|
||||
the interatomic distances are unaffected.
|
||||
</P>
|
||||
<P>The <I>dcd</I> style writes DCD files, a standard atomic trajectory format
|
||||
used by the CHARMM, NAMD, and XPlor molecular dynamics packages. DCD
|
||||
@ -391,7 +402,7 @@ of atom velocity and force and atomic charge.
|
||||
<I>y</I>, <I>z</I> attributes write atom coordinates "unscaled", in the
|
||||
appropriate distance <A HREF = "units.html">units</A> (Angstroms, sigma, etc). Use
|
||||
<I>xs</I>, <I>ys</I>, <I>zs</I> if you want the coordinates "scaled" to the box size,
|
||||
so that each value is 0.0 to 1.0. If the simluation box is triclinic
|
||||
so that each value is 0.0 to 1.0. If the simulation box is triclinic
|
||||
(tilted), then all atom coords will still be between 0.0 and 1.0. Use
|
||||
<I>xu</I>, <I>yu</I>, <I>zu</I> if you want the coordinates "unwrapped" by the image
|
||||
flags for each atom. Unwrapped means that if the atom has passed thru
|
||||
@ -399,7 +410,12 @@ a periodic boundary one or more times, the value is printed for what
|
||||
the coordinate would be if it had not been wrapped back into the
|
||||
periodic box. Note that using <I>xu</I>, <I>yu</I>, <I>zu</I> means that the
|
||||
coordinate values may be far outside the box bounds printed with the
|
||||
snapshot. The image flags can be printed directly using the <I>ix</I>,
|
||||
snapshot. Using <I>xsu</I>, <I>ysu</I>, <I>zsu</I> is similar to using <I>xu</I>, <I>yu</I>, <I>zu</I>,
|
||||
except that the unwrapped coordinates are scaled by the box size. Atoms
|
||||
that have passed through a periodic boundary will have the corresponding
|
||||
cooordinate increased or decreased by 1.0.
|
||||
</P>
|
||||
<P>The image flags can be printed directly using the <I>ix</I>,
|
||||
<I>iy</I>, <I>iz</I> attributes. The <A HREF = "dump_modify.html">dump_modify</A> command
|
||||
describes in more detail what is meant by scaled vs unscaled
|
||||
coordinates and the image flags.
|
||||
|
||||
26
doc/dump.txt
26
doc/dump.txt
@ -37,7 +37,8 @@ args = list of arguments for a particular style :l
|
||||
|
||||
{custom} args = list of atom attributes
|
||||
possible attributes = id, mol, type, mass,
|
||||
x, y, z, xs, ys, zs, xu, yu, zu, ix, iy, iz,
|
||||
x, y, z, xs, ys, zs, xu, yu, zu,
|
||||
xsu, ysu, zsu, ix, iy, iz,
|
||||
vx, vy, vz, fx, fy, fz,
|
||||
q, mux, muy, muz, mu,
|
||||
radius, omegax, omegay, omegaz,
|
||||
@ -52,6 +53,7 @@ args = list of arguments for a particular style :l
|
||||
x,y,z = unscaled atom coordinates
|
||||
xs,ys,zs = scaled atom coordinates
|
||||
xu,yu,zu = unwrapped atom coordinates
|
||||
xsu,ysu,zsu = scaled unwrapped atom coordinates
|
||||
ix,iy,iz = box image that the atom is in
|
||||
vx,vy,vz = atom velocities
|
||||
fx,fy,fz = forces on atoms
|
||||
@ -217,14 +219,23 @@ extended CFG format files, as used by the
|
||||
package. Since the extended CFG format uses a single snapshot of the
|
||||
system per file, a wildcard "*" must be included in the filename, as
|
||||
discussed below. The list of atom attributes for style {cfg} must
|
||||
begin with "id type xs ys zs", since these quantities are needed to
|
||||
begin with either "id type xs ys zs" or "id type xsu ysu zsu" or
|
||||
since these quantities are needed to
|
||||
write the CFG files in the appropriate format (though the "id" and
|
||||
"type" fields do not appear explicitly in the file). Any remaining
|
||||
attributes will be stored as "auxiliary properties" in the CFG files.
|
||||
Note that you will typically want to use the "dump_modify
|
||||
element"_dump_modify.html command with CFG-formatted files, to
|
||||
associate element names with atom types, so that AtomEye can render
|
||||
atoms appropriately.
|
||||
atoms appropriately. When unwrapped coordinates {xsu}, {ysu}, and {zsu}
|
||||
are requested, the nominal AtomEye periodic cell dimensions are expanded
|
||||
by a large factor UNWRAPEXPAND = 10.0, which ensures atoms that are
|
||||
displayed correctly for up to UNWRAPEXPAND/2 periodic boundary crossings
|
||||
in any direction.
|
||||
Beyond this, AtomEye will rewrap the unwrapped coordinates.
|
||||
The expansion causes the atoms to be drawn farther
|
||||
away from the viewer, but it is easy to zoom the atoms closer, and
|
||||
the interatomic distances are unaffected.
|
||||
|
||||
The {dcd} style writes DCD files, a standard atomic trajectory format
|
||||
used by the CHARMM, NAMD, and XPlor molecular dynamics packages. DCD
|
||||
@ -380,7 +391,7 @@ There are several options for outputting atom coordinates. The {x},
|
||||
{y}, {z} attributes write atom coordinates "unscaled", in the
|
||||
appropriate distance "units"_units.html (Angstroms, sigma, etc). Use
|
||||
{xs}, {ys}, {zs} if you want the coordinates "scaled" to the box size,
|
||||
so that each value is 0.0 to 1.0. If the simluation box is triclinic
|
||||
so that each value is 0.0 to 1.0. If the simulation box is triclinic
|
||||
(tilted), then all atom coords will still be between 0.0 and 1.0. Use
|
||||
{xu}, {yu}, {zu} if you want the coordinates "unwrapped" by the image
|
||||
flags for each atom. Unwrapped means that if the atom has passed thru
|
||||
@ -388,7 +399,12 @@ a periodic boundary one or more times, the value is printed for what
|
||||
the coordinate would be if it had not been wrapped back into the
|
||||
periodic box. Note that using {xu}, {yu}, {zu} means that the
|
||||
coordinate values may be far outside the box bounds printed with the
|
||||
snapshot. The image flags can be printed directly using the {ix},
|
||||
snapshot. Using {xsu}, {ysu}, {zsu} is similar to using {xu}, {yu}, {zu},
|
||||
except that the unwrapped coordinates are scaled by the box size. Atoms
|
||||
that have passed through a periodic boundary will have the corresponding
|
||||
cooordinate increased or decreased by 1.0.
|
||||
|
||||
The image flags can be printed directly using the {ix},
|
||||
{iy}, {iz} attributes. The "dump_modify"_dump_modify.html command
|
||||
describes in more detail what is meant by scaled vs unscaled
|
||||
coordinates and the image flags.
|
||||
|
||||
@ -48,14 +48,13 @@ should not be any problems with specifying this fix first in input scripts.
|
||||
<P><I>mode</I> specifies where neighbor list calculations will be performed.
|
||||
If <I>mode</I> is force, neighbor list calculation is performed on the
|
||||
CPU. If <I>mode</I> is force/neigh, neighbor list calculation is
|
||||
performed on the GPU. GPU neighbor
|
||||
list calculation currently cannot be used with a triclinic box.
|
||||
performed on the GPU. GPU neighbor list calculation currently cannot be
|
||||
used with a triclinic box. GPU neighbor list calculation currently
|
||||
cannot be used with <A HREF = "pair_hybrid.html">hybrid</A> pair styles.
|
||||
GPU neighbor lists are not compatible with styles that are not GPU-enabled.
|
||||
When a non-GPU enabled style requires a neighbor list, it will also be
|
||||
built using CPU routines. In these cases, it will typically be more efficient
|
||||
to only use CPU neighbor list builds. For <A HREF = "pair_hybrid.html">hybrid</A> pair
|
||||
styles, GPU calculated neighbor lists might be less efficient because
|
||||
no particles will be skipped in a given neighbor list.
|
||||
to only use CPU neighbor list builds.
|
||||
</P>
|
||||
<P><I>first</I> and <I>last</I> specify the GPUs that will be used for simulation.
|
||||
On each node, the GPU IDs in the inclusive range from <I>first</I> to <I>last</I> will
|
||||
@ -77,7 +76,8 @@ style.
|
||||
</P>
|
||||
<P>In order to use GPU acceleration, a GPU enabled style must be
|
||||
selected in the input script in addition to this fix. Currently,
|
||||
this is limited to a few <A HREF = "pair_style.html">pair styles</A>.
|
||||
this is limited to a few <A HREF = "pair_style.html">pair styles</A> and
|
||||
the PPPM <A HREF = "kspace_style.html">kspace style</A>.
|
||||
</P>
|
||||
<P>More details about these settings and various possible hardware
|
||||
configuration are in <A HREF = "Section_start.html#2_8">this section</A> of the
|
||||
@ -95,8 +95,10 @@ the <A HREF = "run.html">run</A> command.
|
||||
<P><B>Restrictions:</B>
|
||||
</P>
|
||||
<P>The fix must be the first fix specified for a given run. The force/neigh
|
||||
<I>mode</I> should not be used with a triclinic box or GPU-enabled pair styles
|
||||
that need <A HREF = "special_bonds.html">special_bonds</A> settings.
|
||||
<I>mode</I> should not be used with a triclinic box or <A HREF = "pair_hybrid.html">hybrid</A>
|
||||
pair styles.
|
||||
</P>
|
||||
<P><I>split</I> must be positive when using <A HREF = "pair_hybrid.html">hybrid</A> pair styles.
|
||||
</P>
|
||||
<P>Currently, group-ID must be all.
|
||||
</P>
|
||||
|
||||
@ -27,14 +27,21 @@
|
||||
|
||||
<LI>zero or more keyword/value pairs may be appended
|
||||
|
||||
<PRE>keyword = <I>scale</I> or <I>tally</I>
|
||||
<LI>keyword = <I>angmom</I> or <I>omega</I> or <I>scale</I> or <I>tally</I> or <I>zero</I>
|
||||
|
||||
<PRE> <I>angmom</I> value = <I>no</I> or <I>yes</I>
|
||||
<I>no</I> = do not thermostat rotational degrees of freedom via the angular momentum
|
||||
<I>yes</I> = do thermostat rotational degrees of freedom via the angular momentum
|
||||
<I>omega</I> value = <I>no</I> or <I>yes</I>
|
||||
<I>no</I> = do not thermostat rotational degrees of freedom via then angular velocity
|
||||
<I>yes</I> = do thermostat rotational degrees of freedom via the angular velocity
|
||||
<I>scale</I> values = type ratio
|
||||
type = atom type (1-N)
|
||||
ratio = factor by which to scale the damping coefficient
|
||||
<I>tally</I> values = <I>no</I> or <I>yes</I>
|
||||
<I>tally</I> value = <I>no</I> or <I>yes</I>
|
||||
<I>no</I> = do not tally the energy added/subtracted to atoms
|
||||
<I>yes</I> = do tally the energy added/subtracted to atoms
|
||||
<I>zero</I> values = <I>no</I> or <I>yes</I>
|
||||
<I>zero</I> value = <I>no</I> or <I>yes</I>
|
||||
<I>no</I> = do not set total random force to zero
|
||||
<I>yes</I> = set total random force to zero
|
||||
</PRE>
|
||||
@ -135,6 +142,25 @@ generate its own unique seed and its own stream of random numbers.
|
||||
Thus the dynamics of the system will not be identical on two runs on
|
||||
different numbers of processors.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The keyword/value option pairs are used in the following ways.
|
||||
</P>
|
||||
<P>The keyword <I>angmom</I> and <I>omega</I> keywords enable thermostatting of
|
||||
rotational degrees of freedom in addition to the usual translational
|
||||
degrees of freedom. This can only be done for finite-size particles.
|
||||
A simulation using atom_style sphere defines an omega for finite-size
|
||||
spheres. A simulation using atom_style ellipsoid defines a finite
|
||||
size and shape for aspherical particles and an angular momentum. The
|
||||
Langevin formulas for thermostatting the rotational degrees of freedom
|
||||
are the same as those above, where force is replaced by torque, m is
|
||||
replaced by the moment of inertia I, and v is replaced by omega (which
|
||||
is derived from the angular momentum in the case of aspherical
|
||||
particles). The rotational temperature of the particles can be
|
||||
monitored by the <A HREF = "compute_temp_sphere.html">compute temp/sphere</A> and
|
||||
<A HREF = "compute_temp_asphere.html">compute temp/asphere</A> commands with their
|
||||
rotate options.
|
||||
</P>
|
||||
<P>The keyword <I>scale</I> allows the damp factor to be scaled up or down by
|
||||
the specified factor for atoms of that type. This can be useful when
|
||||
different atom types have different sizes or masses. It can be used
|
||||
@ -166,6 +192,8 @@ to zero by subtracting off an equal part of it from each atom in the
|
||||
group. As a result, the center-of-mass of a system with zero initial
|
||||
momentum will not drift over time.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P><B>Restart, fix_modify, output, run start/stop, minimize info:</B>
|
||||
</P>
|
||||
<P>No information about this fix is written to <A HREF = "restart.html">binary restart
|
||||
@ -209,8 +237,8 @@ dpd/tstat</A>
|
||||
</P>
|
||||
<P><B>Default:</B>
|
||||
</P>
|
||||
<P>The option defaults are scale = 1.0 for all types, tally = no, zero =
|
||||
no.
|
||||
<P>The option defaults are angmom = no, omega = no, scale = 1.0 for all
|
||||
types, tally = no, zero = no.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
|
||||
@ -18,14 +18,20 @@ Tstart,Tstop = desired temperature at start/end of run (temperature units) :l
|
||||
damp = damping parameter (time units) :l
|
||||
seed = random number seed to use for white noise (positive integer) :l
|
||||
zero or more keyword/value pairs may be appended :l
|
||||
keyword = {scale} or {tally}
|
||||
keyword = {angmom} or {omega} or {scale} or {tally} or {zero} :l
|
||||
{angmom} value = {no} or {yes}
|
||||
{no} = do not thermostat rotational degrees of freedom via the angular momentum
|
||||
{yes} = do thermostat rotational degrees of freedom via the angular momentum
|
||||
{omega} value = {no} or {yes}
|
||||
{no} = do not thermostat rotational degrees of freedom via then angular velocity
|
||||
{yes} = do thermostat rotational degrees of freedom via the angular velocity
|
||||
{scale} values = type ratio
|
||||
type = atom type (1-N)
|
||||
ratio = factor by which to scale the damping coefficient
|
||||
{tally} values = {no} or {yes}
|
||||
{tally} value = {no} or {yes}
|
||||
{no} = do not tally the energy added/subtracted to atoms
|
||||
{yes} = do tally the energy added/subtracted to atoms
|
||||
{zero} values = {no} or {yes}
|
||||
{zero} value = {no} or {yes}
|
||||
{no} = do not set total random force to zero
|
||||
{yes} = set total random force to zero :pre
|
||||
:ule
|
||||
@ -125,6 +131,25 @@ generate its own unique seed and its own stream of random numbers.
|
||||
Thus the dynamics of the system will not be identical on two runs on
|
||||
different numbers of processors.
|
||||
|
||||
:line
|
||||
|
||||
The keyword/value option pairs are used in the following ways.
|
||||
|
||||
The keyword {angmom} and {omega} keywords enable thermostatting of
|
||||
rotational degrees of freedom in addition to the usual translational
|
||||
degrees of freedom. This can only be done for finite-size particles.
|
||||
A simulation using atom_style sphere defines an omega for finite-size
|
||||
spheres. A simulation using atom_style ellipsoid defines a finite
|
||||
size and shape for aspherical particles and an angular momentum. The
|
||||
Langevin formulas for thermostatting the rotational degrees of freedom
|
||||
are the same as those above, where force is replaced by torque, m is
|
||||
replaced by the moment of inertia I, and v is replaced by omega (which
|
||||
is derived from the angular momentum in the case of aspherical
|
||||
particles). The rotational temperature of the particles can be
|
||||
monitored by the "compute temp/sphere"_compute_temp_sphere.html and
|
||||
"compute temp/asphere"_compute_temp_asphere.html commands with their
|
||||
rotate options.
|
||||
|
||||
The keyword {scale} allows the damp factor to be scaled up or down by
|
||||
the specified factor for atoms of that type. This can be useful when
|
||||
different atom types have different sizes or masses. It can be used
|
||||
@ -156,6 +181,8 @@ to zero by subtracting off an equal part of it from each atom in the
|
||||
group. As a result, the center-of-mass of a system with zero initial
|
||||
momentum will not drift over time.
|
||||
|
||||
:line
|
||||
|
||||
[Restart, fix_modify, output, run start/stop, minimize info:]
|
||||
|
||||
No information about this fix is written to "binary restart
|
||||
@ -199,8 +226,8 @@ dpd/tstat"_pair_dpd.html
|
||||
|
||||
[Default:]
|
||||
|
||||
The option defaults are scale = 1.0 for all types, tally = no, zero =
|
||||
no.
|
||||
The option defaults are angmom = no, omega = no, scale = 1.0 for all
|
||||
types, tally = no, zero = no.
|
||||
|
||||
:line
|
||||
|
||||
|
||||
@ -33,9 +33,13 @@
|
||||
</PRE>
|
||||
<LI>zero or more keyword/value pairs may be appended
|
||||
|
||||
<LI>keyword = <I>temp</I> or <I>press</I> or <I>tparam</I> or <I>pparam</I> or <I>force</I> or <I>torque</I>
|
||||
<LI>keyword = <I>langevin</I> or <I>temp</I> or <I>tparam</I> or <I>force</I> or <I>torque</I>
|
||||
|
||||
<PRE> <I>temp</I> values = Tstart Tstop Tperiod
|
||||
<PRE> <I>langevin</I> values = Tstart Tstop Tperiod seed
|
||||
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
|
||||
Tdamp = temperature damping parameter (time units)
|
||||
seed = random number seed to use for white noise (positive integer)
|
||||
<I>temp</I> values = Tstart Tstop Tdamp
|
||||
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
|
||||
Tdamp = temperature damping parameter (time units)
|
||||
<I>tparam</I> values = Tchain Titer Torder
|
||||
@ -54,7 +58,7 @@
|
||||
<P><B>Examples:</B>
|
||||
</P>
|
||||
<PRE>fix 1 clump rigid single
|
||||
fix 1 clump rigid single force 1 off off on
|
||||
fix 1 clump rigid single force 1 off off on langevin 1.0 1.0 1.0 428984
|
||||
fix 1 polychains rigid/nvt molecule temp 1.0 1.0 5.0
|
||||
fix 1 polychains rigid molecule force 1*5 off off off force 6*10 off off on
|
||||
fix 2 fluid rigid group 3 clump1 clump2 clump3 torque * off off off
|
||||
@ -200,19 +204,35 @@ multiple rigid fixes to be defined, but it is more expensive.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>As stated above, the <I>rigid</I> and <I>rigid/nve</I> styles perform constant
|
||||
NVE time integration. Thus the <I>temp</I>, <I>press</I>, and <I>tparam</I> keywords
|
||||
cannot be used with these styles.
|
||||
<P>The keyword/value option pairs are used in the following ways.
|
||||
</P>
|
||||
<P>The <I>rigid/nvt</I> style performs constant NVT time integration, using a
|
||||
temperature it computes for the rigid bodies which includes their
|
||||
translational and rotational motion. The <I>temp</I> keyword must be used
|
||||
with this style. The desired temperature at each timestep is a ramped
|
||||
value during the run from <I>Tstart</I> to <I>Tstop</I>. The <I>Tdamp</I> parameter
|
||||
is specified in time units and determines how rapidly the temperature
|
||||
is relaxed. For example, a value of 100.0 means to relax the
|
||||
temperature in a timespan of (roughly) 100 time units (tau or fmsec or
|
||||
psec - see the <A HREF = "units.html">units</A> command).
|
||||
<P>The <I>langevin</I> and <I>temp</I> and <I>tparam</I> keywords perform thermostatting
|
||||
of the rigid bodies, altering both their translational and rotational
|
||||
degrees of freedom. What is meant by "temperature" of a collection of
|
||||
rigid bodies and how it can be monitored via the fix output is
|
||||
discussed below.
|
||||
</P>
|
||||
<P>The <I>langevin</I> keyword applies a Langevin thermostat to the constant
|
||||
NVE time integration performed by either the <I>rigid</I> or <I>rigid/nve</I>
|
||||
styles. It cannot be used with the <I>rigid/nvt</I> style. The desired
|
||||
temperature at each timestep is a ramped value during the run from
|
||||
<I>Tstart</I> to <I>Tstop</I>. The <I>Tdamp</I> parameter is specified in time units
|
||||
and determines how rapidly the temperature is relaxed. For example, a
|
||||
value of 100.0 means to relax the temperature in a timespan of
|
||||
(roughly) 100 time units (tau or fmsec or psec - see the
|
||||
<A HREF = "units.html">units</A> command). The random # <I>seed</I> must be a positive
|
||||
integer. The way the Langevin thermostatting operates is explained on
|
||||
the <A HREF = "fix_langevin.html">fix langevin</A> doc page.
|
||||
</P>
|
||||
<P>The <I>temp</I> and <I>tparam</I> keywords apply a Nose/Hoover thermostat to the
|
||||
NVT time integration performed by the <I>rigid/nvt</I> style. They cannot
|
||||
be used with the <I>rigid</I> or <I>rigid/nve</I> styles. The desired
|
||||
temperature at each timestep is a ramped value during the run from
|
||||
<I>Tstart</I> to <I>Tstop</I>. The <I>Tdamp</I> parameter is specified in time units
|
||||
and determines how rapidly the temperature is relaxed. For example, a
|
||||
value of 100.0 means to relax the temperature in a timespan of
|
||||
(roughly) 100 time units (tau or fmsec or psec - see the
|
||||
<A HREF = "units.html">units</A> command).
|
||||
</P>
|
||||
<P>Nose/Hoover chains are used in conjunction with this thermostat. The
|
||||
<I>tparam</I> keyword can optionally be used to change the chain settings
|
||||
@ -222,18 +242,22 @@ oscillations in temperature that can occur in a simulation. As a rule
|
||||
of thumb, increasing the chain length should lead to smaller
|
||||
oscillations.
|
||||
</P>
|
||||
<P>There are alternate ways to thermostat a system of rigid bodies. You
|
||||
can use <A HREF = "fix_langevin.html">fix langevin</A> to treat the system as
|
||||
effectively immersed in an implicit solvent, e.g. a Brownian dynamics
|
||||
model. For hybrid systems with both rigid bodies and solvent
|
||||
particles, you can thermostat only the solvent particles that surround
|
||||
one or more rigid bodies by appropriate choice of groups in the
|
||||
compute and fix commands for temperature and thermostatting. The
|
||||
solvent interactions with the rigid bodies should then effectively
|
||||
thermostat the rigid body temperature as well.
|
||||
<P>IMPORTANT NOTE: There are alternate ways to thermostat a system of
|
||||
rigid bodies. You can use <A HREF = "fix_langevin.html">fix langevin</A> to treat
|
||||
the individual particles in the rigid bodies as effectively immersed
|
||||
in an implicit solvent, e.g. a Brownian dynamics model. For hybrid
|
||||
systems with both rigid bodies and solvent particles, you can
|
||||
thermostat only the solvent particles that surround one or more rigid
|
||||
bodies by appropriate choice of groups in the compute and fix commands
|
||||
for temperature and thermostatting. The solvent interactions with the
|
||||
rigid bodies should then effectively thermostat the rigid body
|
||||
temperature as well without use of the Langevin or Nose/Hoover options
|
||||
associated with the fix rigid commands.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The keyword/value option pairs are used in the following ways.
|
||||
</P>
|
||||
<P>If you use a <A HREF = "compute.html">temperature compute</A> with a group that
|
||||
includes particles in rigid bodies, the degrees-of-freedom removed by
|
||||
each rigid body are accounted for in the temperature (and pressure)
|
||||
@ -289,6 +313,18 @@ rigid/nvt fix to add the energy change induced by the thermostatting
|
||||
to the system's potential energy as part of <A HREF = "thermo_style.html">thermodynamic
|
||||
output</A>.
|
||||
</P>
|
||||
<P>The rigid and rigid/nve fixes computes a global scalar which can be
|
||||
accessed by various <A HREF = "Section_howto.html#4_15">output commands</A>. The
|
||||
scalar value calculated by these fixes is "intensive". The scalar is
|
||||
the current temperature of the collection of rigid bodies. This is
|
||||
averaged over all rigid bodies and their translational and rotational
|
||||
degrees of freedom. The translational energy of a rigid body is 1/2 m
|
||||
v^2, where m = total mass of the body and v = the velocity of its
|
||||
center of mass. The rotational energy of a rigid body is 1/2 I w^2,
|
||||
where I = the moment of inertia tensor of the body and w = its angular
|
||||
velocity. Degrees of freedom constrained by the <I>force</I> and <I>torque</I>
|
||||
keywords are removed from this calculation.
|
||||
</P>
|
||||
<P>The rigid/nvt fix computes a global scalar which can be accessed by
|
||||
various <A HREF = "Section_howto.html#4_15">output commands</A>. The scalar value
|
||||
calculated by the rigid/nvt fix is "extensive". The scalar is the
|
||||
|
||||
@ -24,8 +24,12 @@ bodystyle = {single} or {molecule} or {group} :l
|
||||
groupID1, groupID2, ... = list of N group IDs :pre
|
||||
|
||||
zero or more keyword/value pairs may be appended :l
|
||||
keyword = {temp} or {press} or {tparam} or {pparam} or {force} or {torque} :l
|
||||
{temp} values = Tstart Tstop Tperiod
|
||||
keyword = {langevin} or {temp} or {tparam} or {force} or {torque} :l
|
||||
{langevin} values = Tstart Tstop Tperiod seed
|
||||
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
|
||||
Tdamp = temperature damping parameter (time units)
|
||||
seed = random number seed to use for white noise (positive integer)
|
||||
{temp} values = Tstart Tstop Tdamp
|
||||
Tstart,Tstop = desired temperature at start/stop of run (temperature units)
|
||||
Tdamp = temperature damping parameter (time units)
|
||||
{tparam} values = Tchain Titer Torder
|
||||
@ -43,7 +47,7 @@ keyword = {temp} or {press} or {tparam} or {pparam} or {force} or {torque} :l
|
||||
[Examples:]
|
||||
|
||||
fix 1 clump rigid single
|
||||
fix 1 clump rigid single force 1 off off on
|
||||
fix 1 clump rigid single force 1 off off on langevin 1.0 1.0 1.0 428984
|
||||
fix 1 polychains rigid/nvt molecule temp 1.0 1.0 5.0
|
||||
fix 1 polychains rigid molecule force 1*5 off off off force 6*10 off off on
|
||||
fix 2 fluid rigid group 3 clump1 clump2 clump3 torque * off off off :pre
|
||||
@ -189,19 +193,35 @@ multiple rigid fixes to be defined, but it is more expensive.
|
||||
|
||||
:line
|
||||
|
||||
As stated above, the {rigid} and {rigid/nve} styles perform constant
|
||||
NVE time integration. Thus the {temp}, {press}, and {tparam} keywords
|
||||
cannot be used with these styles.
|
||||
The keyword/value option pairs are used in the following ways.
|
||||
|
||||
The {rigid/nvt} style performs constant NVT time integration, using a
|
||||
temperature it computes for the rigid bodies which includes their
|
||||
translational and rotational motion. The {temp} keyword must be used
|
||||
with this style. The desired temperature at each timestep is a ramped
|
||||
value during the run from {Tstart} to {Tstop}. The {Tdamp} parameter
|
||||
is specified in time units and determines how rapidly the temperature
|
||||
is relaxed. For example, a value of 100.0 means to relax the
|
||||
temperature in a timespan of (roughly) 100 time units (tau or fmsec or
|
||||
psec - see the "units"_units.html command).
|
||||
The {langevin} and {temp} and {tparam} keywords perform thermostatting
|
||||
of the rigid bodies, altering both their translational and rotational
|
||||
degrees of freedom. What is meant by "temperature" of a collection of
|
||||
rigid bodies and how it can be monitored via the fix output is
|
||||
discussed below.
|
||||
|
||||
The {langevin} keyword applies a Langevin thermostat to the constant
|
||||
NVE time integration performed by either the {rigid} or {rigid/nve}
|
||||
styles. It cannot be used with the {rigid/nvt} style. The desired
|
||||
temperature at each timestep is a ramped value during the run from
|
||||
{Tstart} to {Tstop}. The {Tdamp} parameter is specified in time units
|
||||
and determines how rapidly the temperature is relaxed. For example, a
|
||||
value of 100.0 means to relax the temperature in a timespan of
|
||||
(roughly) 100 time units (tau or fmsec or psec - see the
|
||||
"units"_units.html command). The random # {seed} must be a positive
|
||||
integer. The way the Langevin thermostatting operates is explained on
|
||||
the "fix langevin"_fix_langevin.html doc page.
|
||||
|
||||
The {temp} and {tparam} keywords apply a Nose/Hoover thermostat to the
|
||||
NVT time integration performed by the {rigid/nvt} style. They cannot
|
||||
be used with the {rigid} or {rigid/nve} styles. The desired
|
||||
temperature at each timestep is a ramped value during the run from
|
||||
{Tstart} to {Tstop}. The {Tdamp} parameter is specified in time units
|
||||
and determines how rapidly the temperature is relaxed. For example, a
|
||||
value of 100.0 means to relax the temperature in a timespan of
|
||||
(roughly) 100 time units (tau or fmsec or psec - see the
|
||||
"units"_units.html command).
|
||||
|
||||
Nose/Hoover chains are used in conjunction with this thermostat. The
|
||||
{tparam} keyword can optionally be used to change the chain settings
|
||||
@ -211,18 +231,22 @@ oscillations in temperature that can occur in a simulation. As a rule
|
||||
of thumb, increasing the chain length should lead to smaller
|
||||
oscillations.
|
||||
|
||||
There are alternate ways to thermostat a system of rigid bodies. You
|
||||
can use "fix langevin"_fix_langevin.html to treat the system as
|
||||
effectively immersed in an implicit solvent, e.g. a Brownian dynamics
|
||||
model. For hybrid systems with both rigid bodies and solvent
|
||||
particles, you can thermostat only the solvent particles that surround
|
||||
one or more rigid bodies by appropriate choice of groups in the
|
||||
compute and fix commands for temperature and thermostatting. The
|
||||
solvent interactions with the rigid bodies should then effectively
|
||||
thermostat the rigid body temperature as well.
|
||||
IMPORTANT NOTE: There are alternate ways to thermostat a system of
|
||||
rigid bodies. You can use "fix langevin"_fix_langevin.html to treat
|
||||
the individual particles in the rigid bodies as effectively immersed
|
||||
in an implicit solvent, e.g. a Brownian dynamics model. For hybrid
|
||||
systems with both rigid bodies and solvent particles, you can
|
||||
thermostat only the solvent particles that surround one or more rigid
|
||||
bodies by appropriate choice of groups in the compute and fix commands
|
||||
for temperature and thermostatting. The solvent interactions with the
|
||||
rigid bodies should then effectively thermostat the rigid body
|
||||
temperature as well without use of the Langevin or Nose/Hoover options
|
||||
associated with the fix rigid commands.
|
||||
|
||||
:line
|
||||
|
||||
The keyword/value option pairs are used in the following ways.
|
||||
|
||||
If you use a "temperature compute"_compute.html with a group that
|
||||
includes particles in rigid bodies, the degrees-of-freedom removed by
|
||||
each rigid body are accounted for in the temperature (and pressure)
|
||||
@ -278,6 +302,18 @@ rigid/nvt fix to add the energy change induced by the thermostatting
|
||||
to the system's potential energy as part of "thermodynamic
|
||||
output"_thermo_style.html.
|
||||
|
||||
The rigid and rigid/nve fixes computes a global scalar which can be
|
||||
accessed by various "output commands"_Section_howto.html#4_15. The
|
||||
scalar value calculated by these fixes is "intensive". The scalar is
|
||||
the current temperature of the collection of rigid bodies. This is
|
||||
averaged over all rigid bodies and their translational and rotational
|
||||
degrees of freedom. The translational energy of a rigid body is 1/2 m
|
||||
v^2, where m = total mass of the body and v = the velocity of its
|
||||
center of mass. The rotational energy of a rigid body is 1/2 I w^2,
|
||||
where I = the moment of inertia tensor of the body and w = its angular
|
||||
velocity. Degrees of freedom constrained by the {force} and {torque}
|
||||
keywords are removed from this calculation.
|
||||
|
||||
The rigid/nvt fix computes a global scalar which can be accessed by
|
||||
various "output commands"_Section_howto.html#4_15. The scalar value
|
||||
calculated by the rigid/nvt fix is "extensive". The scalar is the
|
||||
|
||||
@ -15,7 +15,7 @@
|
||||
</P>
|
||||
<PRE>kspace_style style value
|
||||
</PRE>
|
||||
<UL><LI>style = <I>none</I> or <I>ewald</I> or <I>pppm</I> or <I>pppm/tip4p</I> or <I>ewald/n</I>
|
||||
<UL><LI>style = <I>none</I> or <I>ewald</I> or <I>pppm</I> or <I>pppm/tip4p</I> or <I>ewald/n</I> or <I>pppm/gpu/single</I> or <I>pppm/gpu/double</I>
|
||||
|
||||
<PRE> <I>none</I> value = none
|
||||
<I>ewald</I> value = precision
|
||||
@ -25,6 +25,10 @@
|
||||
<I>pppm/tip4p</I> value = precision
|
||||
precision = desired accuracy
|
||||
<I>ewald/n</I> value = precision
|
||||
precision = desired accuracy
|
||||
<I>pppm/gpu/single</I> value = precision
|
||||
precision = desired accuracy
|
||||
<I>pppm/gpu/double</I> value = precision
|
||||
precision = desired accuracy
|
||||
</PRE>
|
||||
|
||||
@ -72,6 +76,11 @@ long-range potentials.
|
||||
<P>Currently, only the <I>ewald/n</I> style can be used with non-orthogonal
|
||||
(triclinic symmetry) simulation boxes.
|
||||
</P>
|
||||
<P>The <I>pppm/gpu/single</I> and <I>pppm/gpu/double</I> styles are GPU-enabled
|
||||
version of <I>pppm</I>. See more details below.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>When a kspace style is used, a pair style that includes the
|
||||
short-range correction to the pairwise Coulombic or other 1/r^N forces
|
||||
must also be selected. For Coulombic interactions, these styles are
|
||||
@ -88,6 +97,27 @@ of K-space vectors for style <I>ewald</I> or the FFT grid size for style
|
||||
<P>See the <A HREF = "kspace_modify.html">kspace_modify</A> command for additional
|
||||
options of the K-space solvers that can be set.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The <I>pppm/gpu/single</I> style performs single precision
|
||||
charge assignment and force interpolation calculations on the GPU.
|
||||
The <I>pppm/gpu/double</I> style performs the mesh calculations on the GPU
|
||||
in double precision. FFT solves are calculated on the CPU in both
|
||||
cases. If either <I>pppm/gpu/single</I> or <I>pppm/gpu/double</I> are used with
|
||||
a GPU-enabled pair style, part of the PPPM calculation can be performed
|
||||
concurrently on the GPU while other calculations for non-bonded and
|
||||
bonded force calculation are performed on the CPU.
|
||||
</P>
|
||||
<P>More details about GPU settings and various possible hardware
|
||||
configurations are in <A HREF = "Section_start.html#2_8">this section</A> of the
|
||||
manual.
|
||||
</P>
|
||||
<P>Additional requirements in your input script to run with GPU-enabled
|
||||
PPPM styles are as follows:
|
||||
</P>
|
||||
<P><A HREF = "fix_gpu.html">fix gpu</A> must be used. The fix controls
|
||||
the essential GPU selection and initialization steps.
|
||||
</P>
|
||||
<P><B>Restrictions:</B>
|
||||
</P>
|
||||
<P>A simulation must be 3d and periodic in all dimensions to use an Ewald
|
||||
@ -103,6 +133,11 @@ LAMMPS</A> section for more info.
|
||||
enabled if LAMMPS was built with that package. See the <A HREF = "Section_start.html#2_3">Making
|
||||
LAMMPS</A> section for more info.
|
||||
</P>
|
||||
<P>The <I>pppm/gpu/single</I> and <I>pppm/gpu/double</I> styles are part of the
|
||||
"gpu" package. They are only enabled if LAMMPS was built with that
|
||||
package. See the <A HREF = "Section_start.html#2_3">Making LAMMPS</A> section for
|
||||
more info.
|
||||
</P>
|
||||
<P>When using a long-range pairwise TIP4P potential, you must use kspace
|
||||
style <I>pppm/tip4p</I> and vice versa.
|
||||
</P>
|
||||
|
||||
@ -134,6 +134,7 @@ the pair_style command, and coefficients specified by the associated
|
||||
<LI><A HREF = "pair_lj.html">pair_style lj/cut/coul/long/gpu</A> - GPU-enabled version of LJ with long-range Coulomb
|
||||
<LI><A HREF = "pair_lj.html">pair_style lj/cut/coul/long/tip4p</A> - LJ with long-range Coulomb for TIP4P water
|
||||
<LI><A HREF = "pair_lj_expand.html">pair_style lj/expand</A> - Lennard-Jones for variable size particles
|
||||
<LI><A HREF = "pair_lj_expand.html">pair_style lj/expand/gpu</A> - GPU-enabled version of lj/expand
|
||||
<LI><A HREF = "pair_gromacs.html">pair_style lj/gromacs</A> - GROMACS-style Lennard-Jones potential
|
||||
<LI><A HREF = "pair_gromacs.html">pair_style lj/gromacs/coul/gromacs</A> - GROMACS-style LJ and Coulombic potential
|
||||
<LI><A HREF = "pair_lj_smooth.html">pair_style lj/smooth</A> - smoothed Lennard-Jones potential
|
||||
@ -142,6 +143,7 @@ the pair_style command, and coefficients specified by the associated
|
||||
<LI><A HREF = "pair_lubricate.html">pair_style lubricate</A> - hydrodynamic lubrication forces
|
||||
<LI><A HREF = "pair_meam.html">pair_style meam</A> - modified embedded atom method (MEAM)
|
||||
<LI><A HREF = "pair_morse.html">pair_style morse</A> - Morse potential
|
||||
<LI><A HREF = "pair_morse.html">pair_style morse/gpu</A> - GPU-enabled version of Morse potential
|
||||
<LI><A HREF = "pair_morse.html">pair_style morse/opt</A> - optimized version of Morse potential
|
||||
<LI><A HREF = "pair_peri.html">pair_style peri/lps</A> - peridynamic LPS potential
|
||||
<LI><A HREF = "pair_peri.html">pair_style peri/pmb</A> - peridynamic PMB potential
|
||||
|
||||
@ -11,10 +11,14 @@
|
||||
|
||||
<H3>pair_style lj/expand command
|
||||
</H3>
|
||||
<H3>pair_style lj/expand/gpu command
|
||||
</H3>
|
||||
<P><B>Syntax:</B>
|
||||
</P>
|
||||
<PRE>pair_style lj/expand cutoff
|
||||
</PRE>
|
||||
<PRE>pair_style lj/expand/gpu cutoff
|
||||
</PRE>
|
||||
<UL><LI>cutoff = global cutoff for lj/expand interactions (distance units)
|
||||
</UL>
|
||||
<P><B>Examples:</B>
|
||||
@ -49,6 +53,29 @@ commands, or by mixing as described below:
|
||||
<P>The delta values can be positive or negative. The last coefficient is
|
||||
optional. If not specified, the global LJ cutoff is used.
|
||||
</P>
|
||||
<P>Style <I>lj/expand/gpu</I> is a GPU-enabled version of style <I>lj/expand</I>.
|
||||
See more details below.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The <I>lj/expand/gpu</I> style is identical to the <I>lj/expand</I> style,
|
||||
except that each processor off-loads its pairwise calculations to a
|
||||
GPU chip. Depending on the hardware available on your system this can provide a
|
||||
speed-up. See the <A HREF = "Section_start.html#2_8">Running on GPUs</A> section of
|
||||
the manual for more details about hardware and software requirements
|
||||
for using GPUs.
|
||||
</P>
|
||||
<P>More details about these settings and various possible hardware
|
||||
configuration are in <A HREF = "Section_start.html#2_8">this section</A> of the
|
||||
manual.
|
||||
</P>
|
||||
<P>Additional requirements in your input script to run with GPU-enabled styles
|
||||
are as follows:
|
||||
</P>
|
||||
<P>The <A HREF = "newton.html">newton pair</A> setting must be <I>off</I> and
|
||||
<A HREF = "fix_gpu.html">fix gpu</A> must be used. The fix controls
|
||||
the essential GPU selection and initialization steps.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P><B>Mixing, shift, table, tail correction, restart, rRESPA info</B>:
|
||||
@ -80,7 +107,11 @@ to be specified in an input script that reads a restart file.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P><B>Restrictions:</B> none
|
||||
<P><B>Restrictions:</B>
|
||||
</P>
|
||||
<P>The <I>lj/expand/gpu</I> style is part of the "gpu" package. It is only
|
||||
enabled if LAMMPS was built with that package. See the <A HREF = "Section_start.html#2_3">Making
|
||||
LAMMPS</A> section for more info.
|
||||
</P>
|
||||
<P><B>Related commands:</B>
|
||||
</P>
|
||||
|
||||
@ -11,12 +11,18 @@
|
||||
|
||||
<H3>pair_style morse command
|
||||
</H3>
|
||||
<H3>pair_style morse/gpu command
|
||||
</H3>
|
||||
<H3>pair_style morse/opt command
|
||||
</H3>
|
||||
<P><B>Syntax:</B>
|
||||
</P>
|
||||
<PRE>pair_style morse cutoff
|
||||
</PRE>
|
||||
<PRE>pair_style morse/gpu cutoff
|
||||
</PRE>
|
||||
<PRE>pair_style morse/opt cutoff
|
||||
</PRE>
|
||||
<UL><LI>cutoff = global cutoff for Morse interactions (distance units)
|
||||
</UL>
|
||||
<P><B>Examples:</B>
|
||||
@ -53,6 +59,29 @@ give identical answers. Depending on system size and the processor
|
||||
you are running on, it may be 5-25% faster (for the pairwise portion
|
||||
of the run time).
|
||||
</P>
|
||||
<P>Style <I>morse/gpu</I> is a GPU-enabled version of style <I>morse</I>.
|
||||
See more details below.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P>The <I>morse/gpu</I> style is identical to the <I>morse</I> style,
|
||||
except that each processor off-loads its pairwise calculations to a
|
||||
GPU chip. Depending on the hardware available on your system this can provide a
|
||||
speed-up. See the <A HREF = "Section_start.html#2_8">Running on GPUs</A> section of
|
||||
the manual for more details about hardware and software requirements
|
||||
for using GPUs.
|
||||
</P>
|
||||
<P>More details about these settings and various possible hardware
|
||||
configuration are in <A HREF = "Section_start.html#2_8">this section</A> of the
|
||||
manual.
|
||||
</P>
|
||||
<P>Additional requirements in your input script to run with GPU-enabled styles
|
||||
are as follows:
|
||||
</P>
|
||||
<P>The <A HREF = "newton.html">newton pair</A> setting must be <I>off</I> and
|
||||
<A HREF = "fix_gpu.html">fix gpu</A> must be used. The fix controls
|
||||
the essential GPU selection and initialization steps.
|
||||
</P>
|
||||
<HR>
|
||||
|
||||
<P><B>Mixing, shift, table, tail correction, restart, rRESPA info</B>:
|
||||
@ -82,8 +111,9 @@ to be specified in an input script that reads a restart file.
|
||||
|
||||
<P><B>Restrictions:</B>
|
||||
</P>
|
||||
<P>The <I>morse/opt</I> style is part of the "opt" package. It is only
|
||||
enabled if LAMMPS was built with that package. See the <A HREF = "Section_start.html#2_3">Making
|
||||
<P>The <I>morse/opt</I> style is part of the "opt" package. The <I>morse/gpu</I>
|
||||
style is part of the "gpu" package. They are only
|
||||
enabled if LAMMPS was built with those packages. See the <A HREF = "Section_start.html#2_3">Making
|
||||
LAMMPS</A> section for more info.
|
||||
</P>
|
||||
<P><B>Related commands:</B>
|
||||
|
||||
@ -136,6 +136,7 @@ the pair_style command, and coefficients specified by the associated
|
||||
<LI><A HREF = "pair_lj.html">pair_style lj/cut/coul/long/gpu</A> - GPU-enabled version of LJ with long-range Coulomb
|
||||
<LI><A HREF = "pair_lj.html">pair_style lj/cut/coul/long/tip4p</A> - LJ with long-range Coulomb for TIP4P water
|
||||
<LI><A HREF = "pair_lj_expand.html">pair_style lj/expand</A> - Lennard-Jones for variable size particles
|
||||
<LI><A HREF = "pair_lj_expand.html">pair_style lj/expand/gpu</A> - GPU-enabled version of lj/expand
|
||||
<LI><A HREF = "pair_gromacs.html">pair_style lj/gromacs</A> - GROMACS-style Lennard-Jones potential
|
||||
<LI><A HREF = "pair_gromacs.html">pair_style lj/gromacs/coul/gromacs</A> - GROMACS-style LJ and Coulombic potential
|
||||
<LI><A HREF = "pair_lj_smooth.html">pair_style lj/smooth</A> - smoothed Lennard-Jones potential
|
||||
@ -144,6 +145,7 @@ the pair_style command, and coefficients specified by the associated
|
||||
<LI><A HREF = "pair_lubricate.html">pair_style lubricate</A> - hydrodynamic lubrication forces
|
||||
<LI><A HREF = "pair_meam.html">pair_style meam</A> - modified embedded atom method (MEAM)
|
||||
<LI><A HREF = "pair_morse.html">pair_style morse</A> - Morse potential
|
||||
<LI><A HREF = "pair_morse.html">pair_style morse/gpu</A> - GPU-enabled version of Morse potential
|
||||
<LI><A HREF = "pair_morse.html">pair_style morse/opt</A> - optimized version of Morse potential
|
||||
<LI><A HREF = "pair_peri.html">pair_style peri/lps</A> - peridynamic LPS potential
|
||||
<LI><A HREF = "pair_peri.html">pair_style peri/pmb</A> - peridynamic PMB potential
|
||||
|
||||
@ -1,2 +1,2 @@
|
||||
Geryon Version 10.280
|
||||
|
||||
Geryon Version 11.094
|
||||
|
||||
|
||||
@ -32,13 +32,16 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{ROTATE,ALL};
|
||||
|
||||
#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg != 3 && narg != 4)
|
||||
error->all("Illegal compute temp/asphere command");
|
||||
if (narg < 3) error->all("Illegal compute temp/asphere command");
|
||||
|
||||
scalar_flag = vector_flag = 1;
|
||||
size_vector = 6;
|
||||
@ -48,11 +51,24 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
tempbias = 0;
|
||||
id_bias = NULL;
|
||||
if (narg == 4) {
|
||||
tempbias = 1;
|
||||
int n = strlen(arg[3]) + 1;
|
||||
id_bias = new char[n];
|
||||
strcpy(id_bias,arg[3]);
|
||||
mode = ALL;
|
||||
|
||||
int iarg = 3;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"bias") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal compute temp/asphere command");
|
||||
tempbias = 1;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
id_bias = new char[n];
|
||||
strcpy(id_bias,arg[iarg+1]);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"dof") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal compute temp/asphere command");
|
||||
if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
|
||||
else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
|
||||
else error->all("Illegal compute temp/asphere command");
|
||||
iarg += 2;
|
||||
} else error->all("Illegal compute temp/asphere command");
|
||||
}
|
||||
|
||||
vector = new double[6];
|
||||
@ -76,8 +92,7 @@ ComputeTempAsphere::~ComputeTempAsphere()
|
||||
|
||||
void ComputeTempAsphere::init()
|
||||
{
|
||||
// check that all particles are finite-size
|
||||
// no point particles allowed, spherical is OK
|
||||
// check that all particles are finite-size, no point particles allowed
|
||||
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
@ -114,18 +129,26 @@ void ComputeTempAsphere::init()
|
||||
void ComputeTempAsphere::dof_compute()
|
||||
{
|
||||
// 6 dof for 3d, 3 dof for 2d
|
||||
// which dof are included also depends on mode
|
||||
// assume full rotation of extended particles
|
||||
// user should correct this via compute_modify if needed
|
||||
|
||||
double natoms = group->count(igroup);
|
||||
int nper = 6;
|
||||
if (domain->dimension == 2) nper = 3;
|
||||
int nper;
|
||||
if (domain->dimension == 3) {
|
||||
if (mode == ALL) nper = 6;
|
||||
else nper = 3;
|
||||
} else {
|
||||
if (mode == ALL) nper = 3;
|
||||
else nper = 1;
|
||||
}
|
||||
dof = nper*natoms;
|
||||
|
||||
// additional adjustments to dof
|
||||
|
||||
if (tempbias == 1) dof -= tbias->dof_remove(-1) * natoms;
|
||||
else if (tempbias == 2) {
|
||||
if (tempbias == 1) {
|
||||
if (mode == ALL) dof -= tbias->dof_remove(-1) * natoms;
|
||||
} else if (tempbias == 2) {
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int count = 0;
|
||||
@ -154,46 +177,73 @@ double ComputeTempAsphere::compute_scalar()
|
||||
}
|
||||
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
double **v = atom->v;
|
||||
double **angmom = atom->angmom;
|
||||
double *rmass = atom->rmass;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *shape,*quat;
|
||||
double wbody[3],inertia[3];
|
||||
double rot[3][3];
|
||||
double t = 0.0;
|
||||
|
||||
// sum translationals and rotational energy for each particle
|
||||
// sum translational and rotational energy for each particle
|
||||
// no point particles since divide by inertia
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
double t = 0.0;
|
||||
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
if (mode == ALL) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
|
||||
|
||||
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
|
||||
// principal moments of inertia
|
||||
|
||||
// principal moments of inertia
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
|
||||
inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
|
||||
inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
|
||||
inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
|
||||
inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
t += inertia[0]*wbody[0]*wbody[0] +
|
||||
inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
|
||||
}
|
||||
|
||||
t += inertia[0]*wbody[0]*wbody[0] +
|
||||
inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
// principal moments of inertia
|
||||
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
|
||||
inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
|
||||
inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
|
||||
inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
t += inertia[0]*wbody[0]*wbody[0] +
|
||||
inertia[1]*wbody[1]*wbody[1] + inertia[2]*wbody[2]*wbody[2];
|
||||
}
|
||||
}
|
||||
|
||||
if (tempbias) tbias->restore_bias_all();
|
||||
|
||||
@ -217,58 +267,93 @@ void ComputeTempAsphere::compute_vector()
|
||||
}
|
||||
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
double **v = atom->v;
|
||||
double **angmom = atom->angmom;
|
||||
double *rmass = atom->rmass;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *shape,*quat;
|
||||
double wbody[3],inertia[3];
|
||||
double wbody[3],inertia[3],t[6];
|
||||
double rot[3][3];
|
||||
double massone,t[6];
|
||||
double massone;
|
||||
|
||||
// sum translational and rotational energy for each particle
|
||||
// no point particles since divide by inertia
|
||||
|
||||
for (i = 0; i < 6; i++) t[i] = 0.0;
|
||||
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (mode == ALL) {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
massone = rmass[i];
|
||||
t[0] += massone * v[i][0]*v[i][0];
|
||||
t[1] += massone * v[i][1]*v[i][1];
|
||||
t[2] += massone * v[i][2]*v[i][2];
|
||||
t[3] += massone * v[i][0]*v[i][1];
|
||||
t[4] += massone * v[i][0]*v[i][2];
|
||||
t[5] += massone * v[i][1]*v[i][2];
|
||||
|
||||
// principal moments of inertia
|
||||
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
|
||||
// translational kinetic energy
|
||||
inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
|
||||
inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
|
||||
inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
// rotational kinetic energy
|
||||
|
||||
t[0] += inertia[0]*wbody[0]*wbody[0];
|
||||
t[1] += inertia[1]*wbody[1]*wbody[1];
|
||||
t[2] += inertia[2]*wbody[2]*wbody[2];
|
||||
t[3] += inertia[0]*wbody[0]*wbody[1];
|
||||
t[4] += inertia[1]*wbody[0]*wbody[2];
|
||||
t[5] += inertia[2]*wbody[1]*wbody[2];
|
||||
}
|
||||
|
||||
massone = rmass[i];
|
||||
t[0] += massone * v[i][0]*v[i][0];
|
||||
t[1] += massone * v[i][1]*v[i][1];
|
||||
t[2] += massone * v[i][2]*v[i][2];
|
||||
t[3] += massone * v[i][0]*v[i][1];
|
||||
t[4] += massone * v[i][0]*v[i][2];
|
||||
t[5] += massone * v[i][1]*v[i][2];
|
||||
} else {
|
||||
for (i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
// principal moments of inertia
|
||||
|
||||
// principal moments of inertia
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
massone = rmass[i];
|
||||
|
||||
inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
// rotational kinetic energy
|
||||
|
||||
t[0] += inertia[0]*wbody[0]*wbody[0];
|
||||
t[1] += inertia[1]*wbody[1]*wbody[1];
|
||||
t[2] += inertia[2]*wbody[2]*wbody[2];
|
||||
t[3] += inertia[0]*wbody[0]*wbody[1];
|
||||
t[4] += inertia[1]*wbody[0]*wbody[2];
|
||||
t[5] += inertia[2]*wbody[1]*wbody[2];
|
||||
}
|
||||
inertia[0] = INERTIA*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
|
||||
inertia[1] = INERTIA*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
|
||||
inertia[2] = INERTIA*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat,rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
wbody[0] /= inertia[0];
|
||||
wbody[1] /= inertia[1];
|
||||
wbody[2] /= inertia[2];
|
||||
|
||||
// rotational kinetic energy
|
||||
|
||||
t[0] += inertia[0]*wbody[0]*wbody[0];
|
||||
t[1] += inertia[1]*wbody[1]*wbody[1];
|
||||
t[2] += inertia[2]*wbody[2]*wbody[2];
|
||||
t[3] += inertia[0]*wbody[0]*wbody[1];
|
||||
t[4] += inertia[1]*wbody[0]*wbody[2];
|
||||
t[5] += inertia[2]*wbody[1]*wbody[2];
|
||||
}
|
||||
}
|
||||
|
||||
if (tempbias) tbias->restore_bias_all();
|
||||
|
||||
|
||||
@ -36,7 +36,7 @@ class ComputeTempAsphere : public Compute {
|
||||
void restore_bias(int, double *);
|
||||
|
||||
private:
|
||||
int fix_dof;
|
||||
int fix_dof,mode;
|
||||
double tfactor;
|
||||
char *id_bias;
|
||||
class Compute *tbias; // ptr to additional bias compute
|
||||
|
||||
@ -29,6 +29,8 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
@ -103,9 +105,9 @@ void FixNVEAsphere::initial_integrate(int vflag)
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
|
||||
inertia[0] = rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[1] = rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]) / 5.0;
|
||||
inertia[2] = rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]) / 5.0;
|
||||
inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
|
||||
inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
|
||||
inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
|
||||
|
||||
// compute omega at 1/2 step from angmom at 1/2 step and current q
|
||||
// update quaternion a full step via Richardson iteration
|
||||
|
||||
@ -88,7 +88,9 @@ void FixLangevinEff::post_force_no_tally()
|
||||
f[i][0] += gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
|
||||
f[i][1] += gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
|
||||
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
|
||||
if (abs(spin[i])==1) erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
|
||||
if (abs(spin[i])==1)
|
||||
erforce[i] += 0.75*gamma1*ervel[i] +
|
||||
0.866025404*gamma2*(random->uniform()-0.5);
|
||||
}
|
||||
}
|
||||
} else if (which == BIAS) {
|
||||
@ -105,7 +107,8 @@ void FixLangevinEff::post_force_no_tally()
|
||||
if (v[i][2] != 0.0)
|
||||
f[i][2] += gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
|
||||
if (abs(spin[i])==1 && ervel[i] != 0.0)
|
||||
erforce[i] += 0.75*gamma1*ervel[i] + 0.866025404*gamma2*(random->uniform()-0.5);
|
||||
erforce[i] += 0.75*gamma1*ervel[i] +
|
||||
0.866025404*gamma2*(random->uniform()-0.5);
|
||||
temperature->restore_bias(i,v[i]);
|
||||
}
|
||||
}
|
||||
@ -158,7 +161,8 @@ void FixLangevinEff::post_force_tally()
|
||||
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
|
||||
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
|
||||
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
|
||||
erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
|
||||
erforcelangevin[i] = 0.75*gamma1*ervel[i] +
|
||||
0.866025404*gamma2*(random->uniform()-0.5);
|
||||
f[i][0] += flangevin[i][0];
|
||||
f[i][1] += flangevin[i][1];
|
||||
f[i][2] += flangevin[i][2];
|
||||
@ -175,14 +179,16 @@ void FixLangevinEff::post_force_tally()
|
||||
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
|
||||
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
|
||||
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
|
||||
erforcelangevin[i] = 0.75*gamma1*ervel[i]+0.866025404*gamma2*(random->uniform()-0.5);
|
||||
erforcelangevin[i] = 0.75*gamma1*ervel[i] +
|
||||
0.866025404*gamma2*(random->uniform()-0.5);
|
||||
if (v[i][0] != 0.0) f[i][0] += flangevin[i][0];
|
||||
else flangevin[i][0] = 0.0;
|
||||
if (v[i][1] != 0.0) f[i][1] += flangevin[i][1];
|
||||
else flangevin[i][1] = 0.0;
|
||||
if (v[i][2] != 0.0) f[i][2] += flangevin[i][2];
|
||||
else flangevin[i][2] = 0.0;
|
||||
if (abs(spin[i])==1 && ervel[i] != 0.0) erforce[i] += erforcelangevin[i];
|
||||
if (abs(spin[i])==1 && ervel[i] != 0.0)
|
||||
erforce[i] += erforcelangevin[i];
|
||||
temperature->restore_bias(i,v[i]);
|
||||
}
|
||||
}
|
||||
|
||||
@ -104,7 +104,7 @@ void ComputeClusterAtom::compute_peratom()
|
||||
|
||||
// grow clusterID array if necessary
|
||||
|
||||
if (atom->nlocal > nmax) {
|
||||
if (atom->nlocal+atom->nghost > nmax) {
|
||||
memory->destroy(clusterID);
|
||||
nmax = atom->nmax;
|
||||
memory->create(clusterID,nmax,"cluster/atom:clusterID");
|
||||
|
||||
@ -23,7 +23,7 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define INERTIA 0.4 // moment of inertia for sphere
|
||||
#define INERTIA 0.4 // moment of inertia prefactor for sphere
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
||||
@ -26,15 +26,16 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define INERTIA 0.4 // moment of inertia for sphere
|
||||
enum{ROTATE,ALL};
|
||||
|
||||
#define INERTIA 0.4 // moment of inertia prefactor for sphere
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg)
|
||||
{
|
||||
if (narg != 3 && narg != 4)
|
||||
error->all("Illegal compute temp/sphere command");
|
||||
if (narg < 3) error->all("Illegal compute temp/sphere command");
|
||||
|
||||
scalar_flag = vector_flag = 1;
|
||||
size_vector = 6;
|
||||
@ -44,11 +45,24 @@ ComputeTempSphere::ComputeTempSphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
tempbias = 0;
|
||||
id_bias = NULL;
|
||||
if (narg == 4) {
|
||||
tempbias = 1;
|
||||
int n = strlen(arg[3]) + 1;
|
||||
id_bias = new char[n];
|
||||
strcpy(id_bias,arg[3]);
|
||||
mode = ALL;
|
||||
|
||||
int iarg = 3;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"bias") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal compute temp/sphere command");
|
||||
tempbias = 1;
|
||||
int n = strlen(arg[iarg+1]) + 1;
|
||||
id_bias = new char[n];
|
||||
strcpy(id_bias,arg[iarg+1]);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"dof") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal compute temp/sphere command");
|
||||
if (strcmp(arg[iarg+1],"rotate") == 0) mode = ROTATE;
|
||||
else if (strcmp(arg[iarg+1],"all") == 0) mode = ALL;
|
||||
else error->all("Illegal compute temp/sphere command");
|
||||
iarg += 2;
|
||||
} else error->all("Illegal compute temp/sphere command");
|
||||
}
|
||||
|
||||
vector = new double[6];
|
||||
@ -100,27 +114,34 @@ void ComputeTempSphere::dof_compute()
|
||||
|
||||
// 6 or 3 dof for extended/point particles for 3d
|
||||
// 3 or 2 dof for extended/point particles for 2d
|
||||
// which dof are included also depends on mode
|
||||
// assume full rotation of extended particles
|
||||
// user should correct this via compute_modify if needed
|
||||
|
||||
int dimension = domain->dimension;
|
||||
|
||||
double *radius = atom->radius;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
count = 0;
|
||||
if (dimension == 3) {
|
||||
if (domain->dimension == 3) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (radius[i] == 0.0) count += 3;
|
||||
else count += 6;
|
||||
if (radius[i] == 0.0) {
|
||||
if (mode == ALL) count += 3;
|
||||
} else {
|
||||
if (mode == ALL) count += 6;
|
||||
else count += 3;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (radius[i] == 0.0) count += 2;
|
||||
else count += 3;
|
||||
if (radius[i] == 0.0) {
|
||||
if (mode == ALL) count += 2;
|
||||
} else {
|
||||
if (mode == ALL) count += 3;
|
||||
else count += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -130,28 +151,38 @@ void ComputeTempSphere::dof_compute()
|
||||
// additional adjustments to dof
|
||||
|
||||
if (tempbias == 1) {
|
||||
double natoms = group->count(igroup);
|
||||
dof -= tbias->dof_remove(-1) * natoms;
|
||||
if (mode == ALL) {
|
||||
double natoms = group->count(igroup);
|
||||
dof -= tbias->dof_remove(-1) * natoms;
|
||||
}
|
||||
|
||||
} else if (tempbias == 2) {
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
count = 0;
|
||||
if (dimension == 3) {
|
||||
if (domain->dimension == 3) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (tbias->dof_remove(i)) {
|
||||
if (radius[i] == 0.0) count += 3;
|
||||
else count += 6;
|
||||
if (radius[i] == 0.0) {
|
||||
if (mode == ALL) count += 3;
|
||||
} else {
|
||||
if (mode == ALL) count += 6;
|
||||
else count += 3;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
if (tbias->dof_remove(i)) {
|
||||
if (radius[i] == 0.0) count += 2;
|
||||
else count += 3;
|
||||
if (radius[i] == 0.0) {
|
||||
if (mode == ALL) count += 2;
|
||||
} else {
|
||||
if (mode == ALL) count += 3;
|
||||
else count += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -187,12 +218,19 @@ double ComputeTempSphere::compute_scalar()
|
||||
|
||||
double t = 0.0;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
|
||||
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
|
||||
omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i];
|
||||
}
|
||||
if (mode == ALL) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
|
||||
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
|
||||
omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i];
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
|
||||
omega[i][2]*omega[i][2]) * INERTIA*rmass[i]*radius[i]*radius[i];
|
||||
}
|
||||
|
||||
if (tempbias) tbias->restore_bias_all();
|
||||
|
||||
@ -225,25 +263,38 @@ void ComputeTempSphere::compute_vector()
|
||||
double massone,inertiaone,t[6];
|
||||
for (int i = 0; i < 6; i++) t[i] = 0.0;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
massone = rmass[i];
|
||||
t[0] += massone * v[i][0]*v[i][0];
|
||||
t[1] += massone * v[i][1]*v[i][1];
|
||||
t[2] += massone * v[i][2]*v[i][2];
|
||||
t[3] += massone * v[i][0]*v[i][1];
|
||||
t[4] += massone * v[i][0]*v[i][2];
|
||||
t[5] += massone * v[i][1]*v[i][2];
|
||||
|
||||
inertiaone = INERTIA*radius[i]*radius[i]*rmass[i];
|
||||
t[0] += inertiaone * omega[i][0]*omega[i][0];
|
||||
t[1] += inertiaone * omega[i][1]*omega[i][1];
|
||||
t[2] += inertiaone * omega[i][2]*omega[i][2];
|
||||
t[3] += inertiaone * omega[i][0]*omega[i][1];
|
||||
t[4] += inertiaone * omega[i][0]*omega[i][2];
|
||||
t[5] += inertiaone * omega[i][1]*omega[i][2];
|
||||
}
|
||||
|
||||
if (mode == ALL) {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
massone = rmass[i];
|
||||
t[0] += massone * v[i][0]*v[i][0];
|
||||
t[1] += massone * v[i][1]*v[i][1];
|
||||
t[2] += massone * v[i][2]*v[i][2];
|
||||
t[3] += massone * v[i][0]*v[i][1];
|
||||
t[4] += massone * v[i][0]*v[i][2];
|
||||
t[5] += massone * v[i][1]*v[i][2];
|
||||
|
||||
inertiaone = INERTIA*rmass[i]*radius[i]*radius[i];
|
||||
t[0] += inertiaone * omega[i][0]*omega[i][0];
|
||||
t[1] += inertiaone * omega[i][1]*omega[i][1];
|
||||
t[2] += inertiaone * omega[i][2]*omega[i][2];
|
||||
t[3] += inertiaone * omega[i][0]*omega[i][1];
|
||||
t[4] += inertiaone * omega[i][0]*omega[i][2];
|
||||
t[5] += inertiaone * omega[i][1]*omega[i][2];
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
inertiaone = INERTIA*rmass[i]*radius[i]*radius[i];
|
||||
t[0] += inertiaone * omega[i][0]*omega[i][0];
|
||||
t[1] += inertiaone * omega[i][1]*omega[i][1];
|
||||
t[2] += inertiaone * omega[i][2]*omega[i][2];
|
||||
t[3] += inertiaone * omega[i][0]*omega[i][1];
|
||||
t[4] += inertiaone * omega[i][0]*omega[i][2];
|
||||
t[5] += inertiaone * omega[i][1]*omega[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
if (tempbias) tbias->restore_bias_all();
|
||||
|
||||
MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
@ -36,7 +36,7 @@ class ComputeTempSphere : public Compute {
|
||||
void restore_bias(int, double *);
|
||||
|
||||
private:
|
||||
int fix_dof;
|
||||
int fix_dof,mode;
|
||||
double tfactor;
|
||||
double *inertia;
|
||||
char *id_bias;
|
||||
|
||||
@ -30,6 +30,8 @@
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
#define UNWRAPEXPAND 10.0
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
enum{INT,DOUBLE}; // same as in dump_custom.cpp
|
||||
@ -41,10 +43,20 @@ DumpCFG::DumpCFG(LAMMPS *lmp, int narg, char **arg) :
|
||||
{
|
||||
if (narg < 10 ||
|
||||
strcmp(arg[5],"id") != 0 || strcmp(arg[6],"type") != 0 ||
|
||||
strcmp(arg[7],"xs") != 0 || strcmp(arg[8],"ys") != 0 ||
|
||||
strcmp(arg[9],"zs") != 0)
|
||||
error->all("Dump cfg arguments must start with 'id type xs ys zs'");
|
||||
(strcmp(arg[7],"xs") != 0 && strcmp(arg[7],"xsu") != 0) ||
|
||||
(strcmp(arg[8],"ys") != 0 && strcmp(arg[8],"ysu") != 0) ||
|
||||
(strcmp(arg[9],"zs") != 0 && strcmp(arg[9],"zsu") != 0)
|
||||
)
|
||||
error->all("Dump cfg arguments must start with 'id type xs ys zs' or 'id type xsu ysu zsu'");
|
||||
|
||||
if (strcmp(arg[7],"xs") == 0)
|
||||
if (strcmp(arg[8],"ysu") == 0 || strcmp(arg[9],"zsu") == 0)
|
||||
error->all("Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu");
|
||||
else unwrapflag = 0;
|
||||
else if (strcmp(arg[8],"ys") == 0 || strcmp(arg[9],"zs") == 0)
|
||||
error->all("Dump cfg arguments can not mix xs|ys|zs with xsu|ysu|zsu");
|
||||
else unwrapflag = 1;
|
||||
|
||||
ntypes = atom->ntypes;
|
||||
typenames = NULL;
|
||||
|
||||
@ -189,7 +201,9 @@ void DumpCFG::write_header(bigint n)
|
||||
// special handling for atom style peri
|
||||
// use average volume of particles to scale particles to mimic C atoms
|
||||
// scale box dimension to sc lattice for C with sigma = 1.44 Angstroms
|
||||
|
||||
|
||||
// Special handling for unwrapped coordinates
|
||||
|
||||
double scale;
|
||||
if (atom->peri_flag) {
|
||||
int nlocal = atom->nlocal;
|
||||
@ -199,9 +213,9 @@ void DumpCFG::write_header(bigint n)
|
||||
MPI_Allreduce(&vone,&vave,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
if (atom->natoms) vave /= atom->natoms;
|
||||
if (vave > 0.0) scale = 1.44 / pow(vave,1.0/3.0);
|
||||
else scale = 1.0;
|
||||
} else scale = 1.0;
|
||||
|
||||
} else if (unwrapflag == 1) scale = UNWRAPEXPAND;
|
||||
else scale = 1.0;
|
||||
|
||||
if (me == 0 || multiproc) {
|
||||
char str[64];
|
||||
sprintf(str,"Number of particles = %s\n",BIGINT_FORMAT);
|
||||
@ -261,6 +275,8 @@ void DumpCFG::write_data(int n, double *mybuf)
|
||||
|
||||
// write data lines in rbuf to file after transfer is done
|
||||
|
||||
double unwrap_coord;
|
||||
|
||||
if (nlines == nchosen) {
|
||||
for (itype = 1; itype <= ntypes; itype++) {
|
||||
for (i = 0; i < nchosen; i++)
|
||||
@ -271,11 +287,30 @@ void DumpCFG::write_data(int n, double *mybuf)
|
||||
fprintf(fp,"%s\n",typenames[itype]);
|
||||
for (; i < nchosen; i++) {
|
||||
if (rbuf[i][1] == itype) {
|
||||
for (j = 2; j < size_one; j++) {
|
||||
if (vtype[j] == INT)
|
||||
fprintf(fp,vformat[j],static_cast<int> (rbuf[i][j]));
|
||||
else fprintf(fp,vformat[j],rbuf[i][j]);
|
||||
}
|
||||
if (unwrapflag == 0)
|
||||
for (j = 2; j < size_one; j++) {
|
||||
if (vtype[j] == INT)
|
||||
fprintf(fp,vformat[j],static_cast<int> (rbuf[i][j]));
|
||||
else fprintf(fp,vformat[j],rbuf[i][j]);
|
||||
}
|
||||
else
|
||||
|
||||
// Unwrapped scaled coordinates are shifted to
|
||||
// center of expanded box, to prevent
|
||||
// rewrapping by AtomEye. Dividing by
|
||||
// expansion factor restores correct
|
||||
// interatomic distances.
|
||||
|
||||
for (j = 2; j < 5; j++) {
|
||||
unwrap_coord = (rbuf[i][j] - 0.5)/UNWRAPEXPAND + 0.5;
|
||||
fprintf(fp,vformat[j],unwrap_coord);
|
||||
}
|
||||
for (j = 5; j < size_one; j++) {
|
||||
if (vtype[j] == INT)
|
||||
fprintf(fp,vformat[j],static_cast<int> (rbuf[i][j]));
|
||||
else fprintf(fp,vformat[j],rbuf[i][j]);
|
||||
}
|
||||
|
||||
fprintf(fp,"\n");
|
||||
}
|
||||
}
|
||||
|
||||
@ -36,6 +36,7 @@ class DumpCFG : public DumpCustom {
|
||||
int nchosen; // # of lines to be written on a writing proc
|
||||
int nlines; // # of lines transferred from buf to rbuf
|
||||
double **rbuf; // buf of data lines for data lines rearrangement
|
||||
int unwrapflag; // 1 if unwrapped coordinates are requested
|
||||
|
||||
void init_style();
|
||||
void write_header(bigint);
|
||||
|
||||
@ -35,7 +35,9 @@ using namespace LAMMPS_NS;
|
||||
// same list as in compute_property.cpp, also customize that command
|
||||
|
||||
enum{ID,MOL,TYPE,MASS,
|
||||
X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,IX,IY,IZ,
|
||||
X,Y,Z,XS,YS,ZS,XSTRI,YSTRI,ZSTRI,XU,YU,ZU,XUTRI,YUTRI,ZUTRI,
|
||||
XSU,YSU,ZSU,XSUTRI,YSUTRI,ZSUTRI,
|
||||
IX,IY,IZ,
|
||||
VX,VY,VZ,FX,FY,FZ,
|
||||
Q,MUX,MUY,MUZ,MU,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
|
||||
TQX,TQY,TQZ,SPIN,ERADIUS,ERVEL,ERFORCE,
|
||||
@ -563,6 +565,70 @@ int DumpCustom::count()
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
|
||||
} else if (thresh_array[ithresh] == XSU) {
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
double boxxlo = domain->boxlo[0];
|
||||
double invxprd = 1.0/domain->xprd;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
dchoose[i] = (x[i][0] - boxxlo) * invxprd + (image[i] & 1023) - 512;
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
|
||||
} else if (thresh_array[ithresh] == YSU) {
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
double boxylo = domain->boxlo[1];
|
||||
double invyprd = 1.0/domain->yprd;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
dchoose[i] = (x[i][1] - boxylo) * invyprd + (image[i] >> 10 & 1023) - 512;
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
|
||||
} else if (thresh_array[ithresh] == ZSU) {
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
double boxzlo = domain->boxlo[2];
|
||||
double invzprd = 1.0/domain->zprd;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
dchoose[i] = (x[i][2] - boxzlo) * invzprd + (image[i] >> 20) - 512;
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
|
||||
} else if (thresh_array[ithresh] == XSUTRI) {
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
double *boxlo = domain->boxlo;
|
||||
double *h_inv = domain->h_inv;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
dchoose[i] = h_inv[0]*(x[i][0]-boxlo[0]) +
|
||||
h_inv[5]*(x[i][1]-boxlo[1]) +
|
||||
h_inv[4]*(x[i][2]-boxlo[2]) +
|
||||
(image[i] & 1023) - 512;
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
} else if (thresh_array[ithresh] == YSUTRI) {
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
double *boxlo = domain->boxlo;
|
||||
double *h_inv = domain->h_inv;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
dchoose[i] = h_inv[1]*(x[i][1]-boxlo[1]) +
|
||||
h_inv[3]*(x[i][2]-boxlo[2]) +
|
||||
(image[i] >> 10 & 1023) - 512;
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
} else if (thresh_array[ithresh] == ZSUTRI) {
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
double *boxlo = domain->boxlo;
|
||||
double *h_inv = domain->h_inv;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
dchoose[i] = h_inv[2]*(x[i][2]-boxlo[2]) +
|
||||
(image[i] >> 20) - 512;
|
||||
ptr = dchoose;
|
||||
nstride = 1;
|
||||
|
||||
} else if (thresh_array[ithresh] == IX) {
|
||||
int *image = atom->image;
|
||||
for (i = 0; i < nlocal; i++)
|
||||
@ -879,6 +945,18 @@ void DumpCustom::parse_fields(int narg, char **arg)
|
||||
if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zu_triclinic;
|
||||
else pack_choice[i] = &DumpCustom::pack_zu;
|
||||
vtype[i] = DOUBLE;
|
||||
} else if (strcmp(arg[iarg],"xsu") == 0) {
|
||||
if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_xsu_triclinic;
|
||||
else pack_choice[i] = &DumpCustom::pack_xsu;
|
||||
vtype[i] = DOUBLE;
|
||||
} else if (strcmp(arg[iarg],"ysu") == 0) {
|
||||
if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_ysu_triclinic;
|
||||
else pack_choice[i] = &DumpCustom::pack_ysu;
|
||||
vtype[i] = DOUBLE;
|
||||
} else if (strcmp(arg[iarg],"zsu") == 0) {
|
||||
if (domain->triclinic) pack_choice[i] = &DumpCustom::pack_zsu_triclinic;
|
||||
else pack_choice[i] = &DumpCustom::pack_zsu;
|
||||
vtype[i] = DOUBLE;
|
||||
} else if (strcmp(arg[iarg],"ix") == 0) {
|
||||
pack_choice[i] = &DumpCustom::pack_ix;
|
||||
vtype[i] = INT;
|
||||
@ -1254,6 +1332,19 @@ int DumpCustom::modify_param(int narg, char **arg)
|
||||
else if (strcmp(arg[1],"zu") == 0 && domain->triclinic == 1)
|
||||
thresh_array[nthresh] = ZUTRI;
|
||||
|
||||
else if (strcmp(arg[1],"xsu") == 0 && domain->triclinic == 0)
|
||||
thresh_array[nthresh] = XSU;
|
||||
else if (strcmp(arg[1],"xsu") == 0 && domain->triclinic == 1)
|
||||
thresh_array[nthresh] = XSUTRI;
|
||||
else if (strcmp(arg[1],"ysu") == 0 && domain->triclinic == 0)
|
||||
thresh_array[nthresh] = YSU;
|
||||
else if (strcmp(arg[1],"ysu") == 0 && domain->triclinic == 1)
|
||||
thresh_array[nthresh] = YSUTRI;
|
||||
else if (strcmp(arg[1],"zsu") == 0 && domain->triclinic == 0)
|
||||
thresh_array[nthresh] = ZSU;
|
||||
else if (strcmp(arg[1],"zsu") == 0 && domain->triclinic == 1)
|
||||
thresh_array[nthresh] = ZSUTRI;
|
||||
|
||||
else if (strcmp(arg[1],"ix") == 0) thresh_array[nthresh] = IX;
|
||||
else if (strcmp(arg[1],"iy") == 0) thresh_array[nthresh] = IY;
|
||||
else if (strcmp(arg[1],"iz") == 0) thresh_array[nthresh] = IZ;
|
||||
@ -1821,6 +1912,120 @@ void DumpCustom::pack_zu_triclinic(int n)
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_xsu(int n)
|
||||
{
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double boxxlo = domain->boxlo[0];
|
||||
double invxprd = 1.0/domain->xprd;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (choose[i]) {
|
||||
buf[n] = (x[i][0] - boxxlo) * invxprd + (image[i] & 1023) - 512;
|
||||
n += size_one;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_ysu(int n)
|
||||
{
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double boxylo = domain->boxlo[1];
|
||||
double invyprd = 1.0/domain->yprd;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (choose[i]) {
|
||||
buf[n] = (x[i][1] - boxylo) * invyprd + (image[i] >> 10 & 1023) - 512;
|
||||
n += size_one;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_zsu(int n)
|
||||
{
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double boxzlo = domain->boxlo[2];
|
||||
double invzprd = 1.0/domain->zprd;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (choose[i]) {
|
||||
buf[n] = (x[i][2] - boxzlo) * invzprd + (image[i] >> 20) - 512;
|
||||
n += size_one;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_xsu_triclinic(int n)
|
||||
{
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *boxlo = domain->boxlo;
|
||||
double *h_inv = domain->h_inv;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (choose[i]) {
|
||||
buf[n] = h_inv[0]*(x[i][0]-boxlo[0]) +
|
||||
h_inv[5]*(x[i][1]-boxlo[1]) +
|
||||
h_inv[4]*(x[i][2]-boxlo[2]) +
|
||||
(image[i] & 1023) - 512;
|
||||
n += size_one;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_ysu_triclinic(int n)
|
||||
{
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *boxlo = domain->boxlo;
|
||||
double *h_inv = domain->h_inv;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (choose[i]) {
|
||||
buf[n] = h_inv[1]*(x[i][1]-boxlo[1]) +
|
||||
h_inv[3]*(x[i][2]-boxlo[2]) +
|
||||
(image[i] >> 10 & 1023) - 512;
|
||||
n += size_one;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_zsu_triclinic(int n)
|
||||
{
|
||||
double **x = atom->x;
|
||||
int *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double *boxlo = domain->boxlo;
|
||||
double *h_inv = domain->h_inv;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (choose[i]) {
|
||||
buf[n] = h_inv[2]*(x[i][2]-boxlo[2]) +
|
||||
(image[i] >> 20) - 512;
|
||||
n += size_one;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void DumpCustom::pack_ix(int n)
|
||||
{
|
||||
int *image = atom->image;
|
||||
|
||||
@ -120,6 +120,12 @@ class DumpCustom : public Dump {
|
||||
void pack_xu_triclinic(int);
|
||||
void pack_yu_triclinic(int);
|
||||
void pack_zu_triclinic(int);
|
||||
void pack_xsu(int);
|
||||
void pack_ysu(int);
|
||||
void pack_zsu(int);
|
||||
void pack_xsu_triclinic(int);
|
||||
void pack_ysu_triclinic(int);
|
||||
void pack_zsu_triclinic(int);
|
||||
void pack_ix(int);
|
||||
void pack_iy(int);
|
||||
void pack_iz(int);
|
||||
|
||||
@ -20,7 +20,9 @@
|
||||
#include "string.h"
|
||||
#include "stdlib.h"
|
||||
#include "fix_langevin.h"
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_ellipsoid.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "modify.h"
|
||||
@ -38,6 +40,9 @@ using namespace LAMMPS_NS;
|
||||
|
||||
enum{NOBIAS,BIAS};
|
||||
|
||||
#define SINERTIA 0.4 // moment of inertia prefactor for sphere
|
||||
#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
|
||||
@ -71,12 +76,25 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
|
||||
// optional args
|
||||
|
||||
for (int i = 1; i <= atom->ntypes; i++) ratio[i] = 1.0;
|
||||
oflag = aflag = 0;
|
||||
tally = 0;
|
||||
zeroflag = 0;
|
||||
|
||||
int iarg = 7;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"scale") == 0) {
|
||||
if (strcmp(arg[iarg],"angmom") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal fix langevin command");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) aflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) aflag = 1;
|
||||
else error->all("Illegal fix langevin command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"omega") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal fix langevin command");
|
||||
if (strcmp(arg[iarg+1],"no") == 0) oflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"yes") == 0) oflag = 1;
|
||||
else error->all("Illegal fix langevin command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"scale") == 0) {
|
||||
if (iarg+3 > narg) error->all("Illegal fix langevin command");
|
||||
int itype = atoi(arg[iarg+1]);
|
||||
double scale = atof(arg[iarg+2]);
|
||||
@ -99,6 +117,14 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
|
||||
} else error->all("Illegal fix langevin command");
|
||||
}
|
||||
|
||||
// error check
|
||||
|
||||
if (aflag) {
|
||||
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
||||
if (!avec)
|
||||
error->all("Fix langevin angmom requires atom style ellipsoid");
|
||||
}
|
||||
|
||||
// set temperature = NULL, user can override via fix_modify if wants bias
|
||||
|
||||
id_temp = NULL;
|
||||
@ -140,6 +166,35 @@ int FixLangevin::setmask()
|
||||
|
||||
void FixLangevin::init()
|
||||
{
|
||||
if (oflag && !atom->sphere_flag)
|
||||
error->all("Fix langevin omega require atom style sphere");
|
||||
if (aflag && !atom->ellipsoid_flag)
|
||||
error->all("Fix langevin angmom require atom style ellipsoid");
|
||||
|
||||
// if oflag or aflag set, check that all group particles are finite-size
|
||||
|
||||
if (oflag) {
|
||||
double *radius = atom->radius;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
if (radius[i] == 0.0)
|
||||
error->one("Fix langevin omega requires extended particles");
|
||||
}
|
||||
|
||||
if (aflag) {
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
if (ellipsoid[i] < 0)
|
||||
error->one("Fix langevin angmom requires extended particles");
|
||||
}
|
||||
|
||||
// set force prefactors
|
||||
|
||||
if (!atom->rmass) {
|
||||
@ -219,6 +274,11 @@ void FixLangevin::post_force_no_tally()
|
||||
double fran[3],fsum[3],fsumall[3];
|
||||
fsum[0] = fsum[1] = fsum[2] = 0.0;
|
||||
bigint count;
|
||||
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
|
||||
if (zeroflag) {
|
||||
count = group->count(igroup);
|
||||
@ -227,11 +287,6 @@ void FixLangevin::post_force_no_tally()
|
||||
}
|
||||
|
||||
if (rmass) {
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
|
||||
if (which == NOBIAS) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
@ -280,7 +335,6 @@ void FixLangevin::post_force_no_tally()
|
||||
} else {
|
||||
|
||||
if (which == NOBIAS) {
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
gamma1 = gfactor1[type[i]];
|
||||
@ -295,7 +349,6 @@ void FixLangevin::post_force_no_tally()
|
||||
fsum[1] += fran[1];
|
||||
fsum[2] += fran[2];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
} else if (which == BIAS) {
|
||||
@ -339,6 +392,10 @@ void FixLangevin::post_force_no_tally()
|
||||
}
|
||||
}
|
||||
|
||||
// thermostat omega and angmom
|
||||
|
||||
if (oflag) omega_thermostat(tsqrt);
|
||||
if (aflag) angmom_thermostat(tsqrt);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -374,12 +431,12 @@ void FixLangevin::post_force_tally()
|
||||
// test v = 0 since some computes mask non-participating atoms via v = 0
|
||||
// and added force has extra term not multiplied by v = 0
|
||||
|
||||
if (rmass) {
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
|
||||
if (rmass) {
|
||||
if (which == NOBIAS) {
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
@ -455,6 +512,100 @@ void FixLangevin::post_force_tally()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// thermostat omega and angmom
|
||||
|
||||
if (oflag) omega_thermostat(tsqrt);
|
||||
if (aflag) angmom_thermostat(tsqrt);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
thermostat rotational dof via omega
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixLangevin::omega_thermostat(double tsqrt)
|
||||
{
|
||||
double gamma1,gamma2;
|
||||
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
|
||||
double **torque = atom->torque;
|
||||
double **omega = atom->omega;
|
||||
double *radius = atom->radius;
|
||||
double *rmass = atom->rmass;
|
||||
int *mask = atom->mask;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double tran[3];
|
||||
double inertiaone;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
inertiaone = SINERTIA*radius[i]*radius[i]*rmass[i];
|
||||
gamma1 = -inertiaone / t_period / ftm2v;
|
||||
gamma2 = sqrt(inertiaone) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
|
||||
gamma1 *= 1.0/ratio[type[i]];
|
||||
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
|
||||
tran[0] = gamma2*(random->uniform()-0.5);
|
||||
tran[1] = gamma2*(random->uniform()-0.5);
|
||||
tran[2] = gamma2*(random->uniform()-0.5);
|
||||
torque[i][0] += gamma1*omega[i][0] + tran[0];
|
||||
torque[i][1] += gamma1*omega[i][1] + tran[1];
|
||||
torque[i][2] += gamma1*omega[i][2] + tran[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
thermostat rotational dof via angmom
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixLangevin::angmom_thermostat(double tsqrt)
|
||||
{
|
||||
double gamma1,gamma2;
|
||||
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
double **torque = atom->torque;
|
||||
double **angmom = atom->angmom;
|
||||
double *rmass = atom->rmass;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double inertia[3],wbody[3],omega[3],tran[3],rot[3][3];
|
||||
double *shape,*quat;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
shape = bonus[ellipsoid[i]].shape;
|
||||
inertia[0] = EINERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
|
||||
inertia[1] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
|
||||
inertia[2] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
MathExtra::mq_to_omega(angmom[i],quat,inertia,omega);
|
||||
|
||||
gamma1 = -1.0 / t_period / ftm2v;
|
||||
gamma2 = sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
|
||||
gamma1 *= 1.0/ratio[type[i]];
|
||||
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
|
||||
tran[0] = sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
|
||||
tran[1] = sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
|
||||
tran[2] = sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
|
||||
torque[i][0] += inertia[0]*gamma1*omega[0] + tran[0];
|
||||
torque[i][1] += inertia[1]*gamma1*omega[1] + tran[1];
|
||||
torque[i][2] += inertia[2]*gamma1*omega[2] + tran[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
||||
@ -41,11 +41,13 @@ class FixLangevin : public Fix {
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
int which,tally,zeroflag;
|
||||
int which,tally,zeroflag,oflag,aflag;
|
||||
double t_start,t_stop,t_period;
|
||||
double *gfactor1,*gfactor2,*ratio;
|
||||
double energy,energy_onestep;
|
||||
|
||||
class AtomVecEllipsoid *avec;
|
||||
|
||||
int nmax;
|
||||
double **flangevin;
|
||||
|
||||
@ -57,6 +59,8 @@ class FixLangevin : public Fix {
|
||||
|
||||
virtual void post_force_no_tally();
|
||||
virtual void post_force_tally();
|
||||
void omega_thermostat(double);
|
||||
void angmom_thermostat(double);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -24,7 +24,7 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define INERTIA 0.4 // moment of inertia for sphere
|
||||
#define INERTIA 0.4 // moment of inertia prefactor for sphere
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
||||
@ -24,7 +24,7 @@
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define INERTIA 0.4 // moment of inertia for sphere
|
||||
#define INERTIA 0.4 // moment of inertia prefactor for sphere
|
||||
|
||||
enum{NONE,DIPOLE};
|
||||
|
||||
|
||||
@ -25,6 +25,7 @@
|
||||
#include "modify.h"
|
||||
#include "group.h"
|
||||
#include "comm.h"
|
||||
#include "random_mars.h"
|
||||
#include "force.h"
|
||||
#include "output.h"
|
||||
#include "memory.h"
|
||||
@ -45,11 +46,16 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
{
|
||||
int i,ibody;
|
||||
|
||||
scalar_flag = 1;
|
||||
extscalar = 0;
|
||||
time_integrate = 1;
|
||||
rigid_flag = 1;
|
||||
virial_flag = 1;
|
||||
create_attribute = 1;
|
||||
|
||||
MPI_Comm_rank(world,&me);
|
||||
MPI_Comm_size(world,&nprocs);
|
||||
|
||||
// perform initial allocation of atom-based arrays
|
||||
// register with Atom class
|
||||
|
||||
@ -193,12 +199,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
memory->create(imagebody,nbody,"rigid:imagebody");
|
||||
memory->create(fflag,nbody,3,"rigid:fflag");
|
||||
memory->create(tflag,nbody,3,"rigid:tflag");
|
||||
memory->create(langextra,nbody,6,"rigid:langextra");
|
||||
|
||||
memory->create(sum,nbody,6,"rigid:sum");
|
||||
memory->create(all,nbody,6,"rigid:all");
|
||||
memory->create(remapflag,nbody,4,"rigid:remapflag");
|
||||
|
||||
// initialize force/torque flags to default = 1.0
|
||||
// for 2d: fz, tx, ty = 0.0
|
||||
|
||||
array_flag = 1;
|
||||
size_array_rows = nbody;
|
||||
@ -209,10 +217,13 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
for (i = 0; i < nbody; i++) {
|
||||
fflag[i][0] = fflag[i][1] = fflag[i][2] = 1.0;
|
||||
tflag[i][0] = tflag[i][1] = tflag[i][2] = 1.0;
|
||||
if (domain->dimension == 2) fflag[i][2] = tflag[i][0] = tflag[i][1] = 0.0;
|
||||
}
|
||||
|
||||
// parse optional args
|
||||
|
||||
int seed;
|
||||
langflag = 0;
|
||||
tempflag = 0;
|
||||
pressflag = 0;
|
||||
t_chain = 10;
|
||||
@ -238,6 +249,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
|
||||
else error->all("Illegal fix rigid command");
|
||||
|
||||
if (domain->dimension == 2 && zflag == 1.0)
|
||||
error->all("Fix rigid z force cannot be on for 2d simulation");
|
||||
|
||||
int count = 0;
|
||||
for (int m = mlo; m <= mhi; m++) {
|
||||
fflag[m-1][0] = xflag;
|
||||
@ -266,6 +280,9 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
else if (strcmp(arg[iarg+4],"on") == 0) zflag = 1.0;
|
||||
else error->all("Illegal fix rigid command");
|
||||
|
||||
if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0))
|
||||
error->all("Fix rigid xy torque cannot be on for 2d simulation");
|
||||
|
||||
int count = 0;
|
||||
for (int m = mlo; m <= mhi; m++) {
|
||||
tflag[m-1][0] = xflag;
|
||||
@ -277,10 +294,24 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
iarg += 5;
|
||||
|
||||
} else if (strcmp(arg[iarg],"langevin") == 0) {
|
||||
if (iarg+5 > narg) error->all("Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid") != 0 && strcmp(style,"rigid/nve") != 0)
|
||||
error->all("Illegal fix rigid command");
|
||||
langflag = 1;
|
||||
t_start = atof(arg[iarg+1]);
|
||||
t_stop = atof(arg[iarg+2]);
|
||||
t_period = atof(arg[iarg+3]);
|
||||
seed = atoi(arg[iarg+4]);
|
||||
if (t_period <= 0.0)
|
||||
error->all("Fix rigid langevin period must be > 0.0");
|
||||
if (seed <= 0) error->all("Illegal fix rigid command");
|
||||
iarg += 5;
|
||||
|
||||
} else if (strcmp(arg[iarg],"temp") == 0) {
|
||||
if (iarg+4 > narg) error->all("Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/nvt") != 0 && strcmp(style,"rigid/npt") != 0)
|
||||
error->all("Illegal fix/rigid command");
|
||||
error->all("Illegal fix rigid command");
|
||||
tempflag = 1;
|
||||
t_start = atof(arg[iarg+1]);
|
||||
t_stop = atof(arg[iarg+2]);
|
||||
@ -290,7 +321,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
} else if (strcmp(arg[iarg],"press") == 0) {
|
||||
if (iarg+4 > narg) error->all("Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/npt") != 0)
|
||||
error->all("Illegal fix/rigid command");
|
||||
error->all("Illegal fix rigid command");
|
||||
pressflag = 1;
|
||||
p_start = atof(arg[iarg+1]);
|
||||
p_stop = atof(arg[iarg+2]);
|
||||
@ -300,7 +331,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
} else if (strcmp(arg[iarg],"tparam") == 0) {
|
||||
if (iarg+4 > narg) error->all("Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/nvt") != 0)
|
||||
error->all("Illegal fix/rigid command");
|
||||
error->all("Illegal fix rigid command");
|
||||
t_chain = atoi(arg[iarg+1]);
|
||||
t_iter = atoi(arg[iarg+2]);
|
||||
t_order = atoi(arg[iarg+3]);
|
||||
@ -309,13 +340,18 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
} else if (strcmp(arg[iarg],"pparam") == 0) {
|
||||
if (iarg+2 > narg) error->all("Illegal fix rigid command");
|
||||
if (strcmp(style,"rigid/npt") != 0)
|
||||
error->all("Illegal fix/rigid command");
|
||||
error->all("Illegal fix rigid command");
|
||||
p_chain = atoi(arg[iarg+1]);
|
||||
iarg += 2;
|
||||
|
||||
} else error->all("Illegal fix rigid command");
|
||||
}
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
|
||||
if (langflag) random = new RanMars(lmp,seed + me);
|
||||
else random = NULL;
|
||||
|
||||
// initialize vector output quantities in case accessed before run
|
||||
|
||||
for (i = 0; i < nbody; i++) {
|
||||
@ -369,7 +405,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
|
||||
int nsum = 0;
|
||||
for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody];
|
||||
|
||||
if (comm->me == 0) {
|
||||
if (me == 0) {
|
||||
if (screen) fprintf(screen,"%d rigid bodies with %d atoms\n",nbody,nsum);
|
||||
if (logfile) fprintf(logfile,"%d rigid bodies with %d atoms\n",nbody,nsum);
|
||||
}
|
||||
@ -383,6 +419,8 @@ FixRigid::~FixRigid()
|
||||
|
||||
atom->delete_callback(id,0);
|
||||
|
||||
delete random;
|
||||
|
||||
// delete locally stored arrays
|
||||
|
||||
memory->destroy(body);
|
||||
@ -409,6 +447,7 @@ FixRigid::~FixRigid()
|
||||
memory->destroy(imagebody);
|
||||
memory->destroy(fflag);
|
||||
memory->destroy(tflag);
|
||||
memory->destroy(langextra);
|
||||
|
||||
memory->destroy(sum);
|
||||
memory->destroy(all);
|
||||
@ -422,6 +461,7 @@ int FixRigid::setmask()
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= FINAL_INTEGRATE;
|
||||
if (langflag) mask |= POST_FORCE;
|
||||
mask |= PRE_NEIGHBOR;
|
||||
mask |= INITIAL_INTEGRATE_RESPA;
|
||||
mask |= FINAL_INTEGRATE_RESPA;
|
||||
@ -441,7 +481,7 @@ void FixRigid::init()
|
||||
int count = 0;
|
||||
for (int i = 0; i < modify->nfix; i++)
|
||||
if (strcmp(modify->fix[i]->style,"rigid") == 0) count++;
|
||||
if (count > 1 && comm->me == 0) error->warning("More than one fix rigid");
|
||||
if (count > 1 && me == 0) error->warning("More than one fix rigid");
|
||||
|
||||
// error if npt,nph fix comes before rigid fix
|
||||
|
||||
@ -492,8 +532,8 @@ void FixRigid::init()
|
||||
}
|
||||
|
||||
// grow extended arrays and set extended flags for each particle
|
||||
// dorientflag = 1 if any particles store dipole orientation
|
||||
// qorientflag = 1 if any particles store quat orientation
|
||||
// qorientflag = 1 if any particle stores quat orientation
|
||||
// dorientflag = 1 if any particle stores dipole orientation
|
||||
|
||||
if (extended) {
|
||||
if (atom->mu_flag) dorientflag = 1;
|
||||
@ -639,17 +679,13 @@ void FixRigid::init()
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (body[i] < 0) continue;
|
||||
ibody = body[i];
|
||||
|
||||
itype = type[i];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[itype];
|
||||
massone = rmass[i];
|
||||
|
||||
if (eflags[i] & INERTIA_SPHERE) {
|
||||
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
|
||||
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
|
||||
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
|
||||
}
|
||||
if (eflags[i] & INERTIA_ELLIPSOID) {
|
||||
} else if (eflags[i] & INERTIA_ELLIPSOID) {
|
||||
shape = ebonus[ellipsoid[i]].shape;
|
||||
quatatom = ebonus[ellipsoid[i]].quat;
|
||||
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
|
||||
@ -665,11 +701,12 @@ void FixRigid::init()
|
||||
|
||||
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// diagonalize inertia tensor for each body via Jacobi rotations
|
||||
// inertia = 3 eigenvalues = principal moments of inertia
|
||||
// ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
|
||||
|
||||
// evectors and exzy_space = 3 evectors = principal axes of rigid body
|
||||
|
||||
int ierror;
|
||||
double ez0,ez1,ez2;
|
||||
double cross[3];
|
||||
double tensor[3][3],evectors[3][3];
|
||||
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
@ -686,11 +723,9 @@ void FixRigid::init()
|
||||
ex_space[ibody][0] = evectors[0][0];
|
||||
ex_space[ibody][1] = evectors[1][0];
|
||||
ex_space[ibody][2] = evectors[2][0];
|
||||
|
||||
ey_space[ibody][0] = evectors[0][1];
|
||||
ey_space[ibody][1] = evectors[1][1];
|
||||
ey_space[ibody][2] = evectors[2][1];
|
||||
|
||||
ez_space[ibody][0] = evectors[0][2];
|
||||
ez_space[ibody][1] = evectors[1][2];
|
||||
ez_space[ibody][2] = evectors[2][2];
|
||||
@ -706,21 +741,11 @@ void FixRigid::init()
|
||||
if (inertia[ibody][2] < EPSILON*max) inertia[ibody][2] = 0.0;
|
||||
|
||||
// enforce 3 evectors as a right-handed coordinate system
|
||||
// flip 3rd evector if needed
|
||||
|
||||
ez0 = ex_space[ibody][1]*ey_space[ibody][2] -
|
||||
ex_space[ibody][2]*ey_space[ibody][1];
|
||||
ez1 = ex_space[ibody][2]*ey_space[ibody][0] -
|
||||
ex_space[ibody][0]*ey_space[ibody][2];
|
||||
ez2 = ex_space[ibody][0]*ey_space[ibody][1] -
|
||||
ex_space[ibody][1]*ey_space[ibody][0];
|
||||
|
||||
if (ez0*ez_space[ibody][0] + ez1*ez_space[ibody][1] +
|
||||
ez2*ez_space[ibody][2] < 0.0) {
|
||||
ez_space[ibody][0] = -ez_space[ibody][0];
|
||||
ez_space[ibody][1] = -ez_space[ibody][1];
|
||||
ez_space[ibody][2] = -ez_space[ibody][2];
|
||||
}
|
||||
// flip 3rd vector if needed
|
||||
|
||||
MathExtra::cross3(ex_space[ibody],ey_space[ibody],cross);
|
||||
if (MathExtra::dot3(cross,ez_space[ibody]) < 0.0)
|
||||
MathExtra::negate3(ez_space[ibody]);
|
||||
|
||||
// create initial quaternion
|
||||
|
||||
@ -814,17 +839,13 @@ void FixRigid::init()
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (body[i] < 0) continue;
|
||||
ibody = body[i];
|
||||
|
||||
itype = type[i];
|
||||
if (rmass) massone = rmass[i];
|
||||
else massone = mass[itype];
|
||||
massone = rmass[i];
|
||||
|
||||
if (eflags[i] & INERTIA_SPHERE) {
|
||||
sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
|
||||
sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
|
||||
sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
|
||||
}
|
||||
if (eflags[i] & INERTIA_ELLIPSOID) {
|
||||
} else if (eflags[i] & INERTIA_ELLIPSOID) {
|
||||
shape = ebonus[ellipsoid[i]].shape;
|
||||
MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
|
||||
sum[ibody][0] += ivec[0];
|
||||
@ -868,6 +889,15 @@ void FixRigid::init()
|
||||
fabs(all[ibody][5]/norm) > TOLERANCE)
|
||||
error->all("Fix rigid: Bad principal moments");
|
||||
}
|
||||
|
||||
// temperature scale factor
|
||||
|
||||
double ndof = 0.0;
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
|
||||
ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
|
||||
}
|
||||
tfactor = force->mvv2e / (ndof * force->boltz);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -1011,6 +1041,13 @@ void FixRigid::setup(int vflag)
|
||||
torque[ibody][2] = all[ibody][5];
|
||||
}
|
||||
|
||||
// zero langextra in case Langevin thermostat not used
|
||||
// no point to calling post_force() here since langextra
|
||||
// is only added to fcm/torque in final_integrate()
|
||||
|
||||
for (ibody = 0; ibody < nbody; ibody++)
|
||||
for (i = 0; i < 6; i++) langextra[ibody][i] = 0.0;
|
||||
|
||||
// virial setup before call to set_v
|
||||
|
||||
if (vflag) v_setup(vflag);
|
||||
@ -1085,6 +1122,50 @@ void FixRigid::initial_integrate(int vflag)
|
||||
set_xv();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply Langevin thermostat to all 6 DOF of rigid bodies
|
||||
computed by proc 0, broadcast to other procs
|
||||
unlike fix langevin, this stores extra force in extra arrays,
|
||||
which are added in when final_integrate() calculates a new fcm/torque
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixRigid::post_force(int vflag)
|
||||
{
|
||||
if (me == 0) {
|
||||
double gamma1,gamma2;
|
||||
|
||||
double delta = update->ntimestep - update->beginstep;
|
||||
delta /= update->endstep - update->beginstep;
|
||||
double t_target = t_start + delta * (t_stop-t_start);
|
||||
double tsqrt = sqrt(t_target);
|
||||
|
||||
double boltz = force->boltz;
|
||||
double dt = update->dt;
|
||||
double mvv2e = force->mvv2e;
|
||||
double ftm2v = force->ftm2v;
|
||||
|
||||
for (int i = 0; i < nbody; i++) {
|
||||
gamma1 = -masstotal[i] / t_period / ftm2v;
|
||||
gamma2 = sqrt(masstotal[i]) * tsqrt *
|
||||
sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
|
||||
langextra[i][0] = gamma1*vcm[i][0] + gamma2*(random->uniform()-0.5);
|
||||
langextra[i][1] = gamma1*vcm[i][1] + gamma2*(random->uniform()-0.5);
|
||||
langextra[i][2] = gamma1*vcm[i][2] + gamma2*(random->uniform()-0.5);
|
||||
|
||||
gamma1 = -1.0 / t_period / ftm2v;
|
||||
gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
|
||||
langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] +
|
||||
sqrt(inertia[i][0])*gamma2*(random->uniform()-0.5);
|
||||
langextra[i][4] = inertia[i][1]*gamma1*omega[i][1] +
|
||||
sqrt(inertia[i][1])*gamma2*(random->uniform()-0.5);
|
||||
langextra[i][5] = inertia[i][2]*gamma1*omega[i][2] +
|
||||
sqrt(inertia[i][2])*gamma2*(random->uniform()-0.5);
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Bcast(&langextra[0][0],6*nbody,MPI_DOUBLE,0,world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixRigid::final_integrate()
|
||||
@ -1163,13 +1244,17 @@ void FixRigid::final_integrate()
|
||||
|
||||
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// update vcm and angmom
|
||||
// include Langevin thermostat forces
|
||||
// fflag,tflag = 0 for some dimensions in 2d
|
||||
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
fcm[ibody][0] = all[ibody][0];
|
||||
fcm[ibody][1] = all[ibody][1];
|
||||
fcm[ibody][2] = all[ibody][2];
|
||||
torque[ibody][0] = all[ibody][3];
|
||||
torque[ibody][1] = all[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5];
|
||||
fcm[ibody][0] = all[ibody][0] + langextra[ibody][0];
|
||||
fcm[ibody][1] = all[ibody][1] + langextra[ibody][1];
|
||||
fcm[ibody][2] = all[ibody][2] + langextra[ibody][2];
|
||||
torque[ibody][0] = all[ibody][3] + langextra[ibody][3];
|
||||
torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
|
||||
|
||||
// update vcm by 1/2 step
|
||||
|
||||
@ -1373,7 +1458,7 @@ int FixRigid::dof(int igroup)
|
||||
if (nall[ibody]+mall[ibody] > 0 &&
|
||||
nall[ibody]+mall[ibody] != nrigid[ibody]) flag = 1;
|
||||
}
|
||||
if (flag && comm->me == 0)
|
||||
if (flag && me == 0)
|
||||
error->warning("Computing temperature of portions of rigid bodies");
|
||||
|
||||
// remove appropriate DOFs for each rigid body wholly in temperature group
|
||||
@ -1847,6 +1932,42 @@ void FixRigid::reset_dt()
|
||||
dtq = 0.5 * update->dt;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return temperature of collection of rigid bodies
|
||||
non-active DOF are removed by fflag/tflag and in tfactor
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double FixRigid::compute_scalar()
|
||||
{
|
||||
double wbody[3],rot[3][3];
|
||||
|
||||
double t = 0.0;
|
||||
|
||||
for (int i = 0; i < nbody; i++) {
|
||||
t += masstotal[i] * (fflag[i][0]*vcm[i][0]*vcm[i][0] +
|
||||
fflag[i][1]*vcm[i][1]*vcm[i][1] + \
|
||||
fflag[i][2]*vcm[i][2]*vcm[i][2]);
|
||||
|
||||
// wbody = angular velocity in body frame
|
||||
|
||||
MathExtra::quat_to_mat(quat[i],rot);
|
||||
MathExtra::transpose_matvec(rot,angmom[i],wbody);
|
||||
if (inertia[i][0] == 0.0) wbody[0] = 0.0;
|
||||
else wbody[0] /= inertia[i][0];
|
||||
if (inertia[i][1] == 0.0) wbody[1] = 0.0;
|
||||
else wbody[1] /= inertia[i][1];
|
||||
if (inertia[i][2] == 0.0) wbody[2] = 0.0;
|
||||
else wbody[2] /= inertia[i][2];
|
||||
|
||||
t += tflag[i][0]*inertia[i][0]*wbody[0]*wbody[0] +
|
||||
tflag[i][1]*inertia[i][1]*wbody[1]*wbody[1] +
|
||||
tflag[i][2]*inertia[i][2]*wbody[2]*wbody[2];
|
||||
}
|
||||
|
||||
t *= tfactor;
|
||||
return t;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
return attributes of a rigid body
|
||||
15 values per body
|
||||
|
||||
@ -32,9 +32,11 @@ class FixRigid : public Fix {
|
||||
virtual void init();
|
||||
virtual void setup(int);
|
||||
virtual void initial_integrate(int);
|
||||
void post_force(int);
|
||||
virtual void final_integrate();
|
||||
void initial_integrate_respa(int, int, int);
|
||||
void final_integrate_respa(int, int);
|
||||
virtual double compute_scalar();
|
||||
|
||||
double memory_usage();
|
||||
void grow_arrays(int);
|
||||
@ -50,6 +52,7 @@ class FixRigid : public Fix {
|
||||
double compute_array(int, int);
|
||||
|
||||
protected:
|
||||
int me,nprocs;
|
||||
double dtv,dtf,dtq;
|
||||
double *step_respa;
|
||||
int triclinic;
|
||||
@ -70,6 +73,7 @@ class FixRigid : public Fix {
|
||||
int *imagebody; // image flags of xcm of each rigid body
|
||||
double **fflag; // flag for on/off of center-of-mass force
|
||||
double **tflag; // flag for on/off of center-of-mass torque
|
||||
double **langextra; // Langevin thermostat forces and torques
|
||||
|
||||
int *body; // which body each atom is part of (-1 if none)
|
||||
double **displace; // displacement of each atom in body coords
|
||||
@ -85,6 +89,9 @@ class FixRigid : public Fix {
|
||||
double **qorient; // rotation state of ext particle wrt rigid body
|
||||
double **dorient; // orientation of dipole mu wrt rigid body
|
||||
|
||||
double tfactor; // scale factor on temperature of rigid bodies
|
||||
int langflag; // 0/1 = no/yes Langevin thermostat
|
||||
|
||||
int tempflag; // NVT settings
|
||||
double t_start,t_stop;
|
||||
double t_period,t_freq;
|
||||
@ -95,6 +102,7 @@ class FixRigid : public Fix {
|
||||
double p_period,p_freq;
|
||||
int p_chain;
|
||||
|
||||
class RanMars *random;
|
||||
class AtomVecEllipsoid *avec_ellipsoid;
|
||||
|
||||
// bitmasks for eflags
|
||||
|
||||
@ -223,16 +223,20 @@ void FixRigidNVE::final_integrate()
|
||||
|
||||
MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
|
||||
|
||||
// update vcm and angmom
|
||||
// include Langevin thermostat forces
|
||||
// fflag,tflag = 0 for some dimensions in 2d
|
||||
|
||||
double mbody[3],tbody[3],fquat[4];
|
||||
double dtf2 = dtf * 2.0;
|
||||
|
||||
|
||||
for (ibody = 0; ibody < nbody; ibody++) {
|
||||
fcm[ibody][0] = all[ibody][0];
|
||||
fcm[ibody][1] = all[ibody][1];
|
||||
fcm[ibody][2] = all[ibody][2];
|
||||
torque[ibody][0] = all[ibody][3];
|
||||
torque[ibody][1] = all[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5];
|
||||
fcm[ibody][0] = all[ibody][0] + langextra[ibody][0];
|
||||
fcm[ibody][1] = all[ibody][1] + langextra[ibody][1];
|
||||
fcm[ibody][2] = all[ibody][2] + langextra[ibody][2];
|
||||
torque[ibody][0] = all[ibody][3] + langextra[ibody][3];
|
||||
torque[ibody][1] = all[ibody][4] + langextra[ibody][4];
|
||||
torque[ibody][2] = all[ibody][5] + langextra[ibody][5];
|
||||
|
||||
// update vcm by 1/2 step
|
||||
|
||||
|
||||
@ -177,7 +177,7 @@ void rotate(double matrix[3][3], int i, int j, int k, int l,
|
||||
/* ----------------------------------------------------------------------
|
||||
Richardson iteration to update quaternion from angular momentum
|
||||
return new normalized quaternion q
|
||||
also returns
|
||||
also returns updated omega at 1/2 step
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void richardson(double *q, double *m, double *w, double *moments, double dtq)
|
||||
@ -487,7 +487,7 @@ void inertia_line(double length, double theta, double mass, double *inertia)
|
||||
/* ----------------------------------------------------------------------
|
||||
compute space-frame inertia tensor of a triangle
|
||||
v0,v1,v2 = 3 vertices of triangle
|
||||
from http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle:
|
||||
from http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle
|
||||
inertia tensor = a/24 (v0^2 + v1^2 + v2^2 + (v0+v1+v2)^2) I - a Vt S V
|
||||
a = 2*area of tri = |(v1-v0) x (v2-v0)|
|
||||
I = 3x3 identity matrix
|
||||
@ -506,9 +506,9 @@ void inertia_triangle(double *v0, double *v1, double *v2,
|
||||
double v[3][3],sv[3][3],vtsv[3][3];
|
||||
double vvv[3],v1mv0[3],v2mv0[3],normal[3];
|
||||
|
||||
v[0][0] = v0[0]; v[0][1] = v0[2]; v[0][2] = v0[3];
|
||||
v[1][0] = v1[0]; v[1][1] = v1[2]; v[1][2] = v1[3];
|
||||
v[2][0] = v2[0]; v[2][1] = v2[2]; v[2][2] = v2[3];
|
||||
v[0][0] = v0[0]; v[0][1] = v0[1]; v[0][2] = v0[2];
|
||||
v[1][0] = v1[0]; v[1][1] = v1[1]; v[1][2] = v1[2];
|
||||
v[2][0] = v2[0]; v[2][1] = v2[1]; v[2][2] = v2[2];
|
||||
|
||||
times3(s,v,sv);
|
||||
transpose_times3(v,sv,vtsv);
|
||||
@ -523,9 +523,7 @@ void inertia_triangle(double *v0, double *v1, double *v2,
|
||||
sub3(v2,v0,v2mv0);
|
||||
cross3(v1mv0,v2mv0,normal);
|
||||
double a = len3(normal);
|
||||
double inv24 = 1.0/24.0;
|
||||
|
||||
// NOTE: use mass
|
||||
double inv24 = mass/24.0;
|
||||
|
||||
inertia[0] = inv24*a*(sum-vtsv[0][0]);
|
||||
inertia[1] = inv24*a*(sum-vtsv[1][1]);
|
||||
@ -535,6 +533,30 @@ void inertia_triangle(double *v0, double *v1, double *v2,
|
||||
inertia[5] = -inv24*a*vtsv[0][1];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
compute space-frame inertia tensor of a triangle
|
||||
idiag = previously computed diagonal inertia tensor
|
||||
quat = orientiation quaternion of triangle
|
||||
return symmetric inertia tensor as 6-vector in Voigt notation
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void inertia_triangle(double *idiag, double *quat, double mass,
|
||||
double *inertia)
|
||||
{
|
||||
double p[3][3],ptrans[3][3],itemp[3][3],tensor[3][3];
|
||||
|
||||
quat_to_mat(quat,p);
|
||||
quat_to_mat_trans(quat,ptrans);
|
||||
diag_times3(idiag,ptrans,itemp);
|
||||
times3(p,itemp,tensor);
|
||||
inertia[0] = tensor[0][0];
|
||||
inertia[1] = tensor[1][1];
|
||||
inertia[2] = tensor[2][2];
|
||||
inertia[3] = tensor[1][2];
|
||||
inertia[4] = tensor[0][2];
|
||||
inertia[5] = tensor[0][1];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
}
|
||||
|
||||
@ -30,6 +30,8 @@ namespace MathExtra {
|
||||
inline void norm3(double *v);
|
||||
inline void normalize3(const double *v, double *ans);
|
||||
inline void snormalize3(const double, const double *v, double *ans);
|
||||
inline void negate3(double *v);
|
||||
inline void scale3(double s, double *v);
|
||||
inline void add3(const double *v1, const double *v2, double *ans);
|
||||
inline void sub3(const double *v1, const double *v2, double *ans);
|
||||
inline double len3(const double *v);
|
||||
@ -118,6 +120,8 @@ namespace MathExtra {
|
||||
double *inertia);
|
||||
void inertia_triangle(double *v0, double *v1, double *v2,
|
||||
double mass, double *inertia);
|
||||
void inertia_triangle(double *idiag, double *quat, double mass,
|
||||
double *inertia);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -156,6 +160,28 @@ void MathExtra::snormalize3(const double length, const double *v, double *ans)
|
||||
ans[2] = v[2]*scale;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
negate vector v
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MathExtra::negate3(double *v)
|
||||
{
|
||||
v[0] = -v[0];
|
||||
v[1] = -v[1];
|
||||
v[2] = -v[2];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
scale vector v by s
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void MathExtra::scale3(double s, double *v)
|
||||
{
|
||||
v[0] *= s;
|
||||
v[1] *= s;
|
||||
v[2] *= s;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
ans = v1 + v2
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
@ -45,7 +45,7 @@ class Memory : protected Pointers {
|
||||
bigint nbytes = sizeof(TYPE) * n;
|
||||
array = (TYPE *) smalloc(nbytes,name);
|
||||
return array;
|
||||
};
|
||||
}
|
||||
|
||||
template <typename TYPE>
|
||||
TYPE **create(TYPE **&array, int n, const char *name) {fail(name);}
|
||||
@ -62,7 +62,7 @@ class Memory : protected Pointers {
|
||||
bigint nbytes = sizeof(TYPE) * n;
|
||||
array = (TYPE *) srealloc(array,nbytes,name);
|
||||
return array;
|
||||
};
|
||||
}
|
||||
|
||||
template <typename TYPE>
|
||||
TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);}
|
||||
@ -75,7 +75,7 @@ class Memory : protected Pointers {
|
||||
void destroy(TYPE *array)
|
||||
{
|
||||
sfree(array);
|
||||
};
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
create a 1d array with index from nlo to nhi inclusive
|
||||
|
||||
@ -1 +1 @@
|
||||
#define LAMMPS_VERSION "20 Apr 2011"
|
||||
#define LAMMPS_VERSION "2 May 2011"
|
||||
|
||||
Reference in New Issue
Block a user