Merge branch 'master' into lammps-icms

This commit is contained in:
Axel Kohlmeyer
2010-11-25 10:24:09 -05:00
39 changed files with 603 additions and 177 deletions

Binary file not shown.

Before

Width:  |  Height:  |  Size: 27 KiB

After

Width:  |  Height:  |  Size: 7.3 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 20 KiB

After

Width:  |  Height:  |  Size: 5.7 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 27 KiB

After

Width:  |  Height:  |  Size: 7.5 KiB

BIN
doc/JPG/screenshot_vmd.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 320 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.6 KiB

View File

@ -609,14 +609,14 @@ distribution.
<TR><TD >gui.py</TD><TD > GUI go/stop/temperature-slider to control LAMMPS</TD></TR>
<TR><TD >plot.py</TD><TD > real-time temeperature plot with GnuPlot via Pizza.py</TD></TR>
<TR><TD >viz_tool.py</TD><TD > real-time viz via some viz package</TD></TR>
<TR><TD >vizplotgui_tool.py</TD><TD > combination of viz.py and plot.py and gui.py
<TR><TD >vizplotgui_tool.py</TD><TD > combination of viz_tool.py and plot.py and gui.py
</TD></TR></TABLE></DIV>
<HR>
<P>For the viz_tool.py and vizplotgui_tool.py commands, replace "tool"
with "gl" or "atomeye" or "pymol", depending on what visualization
package you have installed. We hope to add a VMD option soon.
with "gl" or "atomeye" or "pymol" or "vmd", depending on what
visualization package you have installed.
</P>
<P>Note that for GL, you need to be able to run the Pizza.py GL tool,
which is included in the pizza sub-directory. See the <A HREF = "http://www.sandia.gov/~sjplimp/pizza.html">Pizza.py doc
@ -624,7 +624,7 @@ pages</A> for more info:
</P>
<P>Note that for AtomEye, you need version 3, and their is a line in the
<P>Note that for AtomEye, you need version 3, and there is a line in the
scripts that specifies the path and name of the executable. See the
AtomEye WWW pages <A HREF = "http://mt.seas.upenn.edu/Archive/Graphics/A">here</A> or <A HREF = "http://mt.seas.upenn.edu/Archive/Graphics/A3/A3.html">here</A> for more details:
</P>
@ -652,6 +652,10 @@ http://sourceforge.net/scm/?type=svn&group_id=4546
<P>The latter link is to the open-source version.
</P>
<P>Note that for VMD, you need a fairly current version (1.8.7 works for
me) and there are some lines in the pizza/vmd.py script for 4 PIZZA
variables that have to match the VMD installation on your system.
</P>
<HR>
<P>See the python/README file for instructions on how to run them and the
@ -666,4 +670,6 @@ different visualization package options. Click to see larger images:
<A HREF = "JPG/screenshot_pymol.jpg"><IMG SRC = "JPG/screenshot_pymol_small.jpg"></A>
<A HREF = "JPG/screenshot_vmd.jpg"><IMG SRC = "JPG/screenshot_vmd_small.jpg"></A>
</HTML>

View File

@ -602,13 +602,13 @@ simple.py, mimic operation of couple/simple/simple.cpp in Python,
gui.py, GUI go/stop/temperature-slider to control LAMMPS,
plot.py, real-time temeperature plot with GnuPlot via Pizza.py,
viz_tool.py, real-time viz via some viz package,
vizplotgui_tool.py, combination of viz.py and plot.py and gui.py :tb(c=2)
vizplotgui_tool.py, combination of viz_tool.py and plot.py and gui.py :tb(c=2)
:line
For the viz_tool.py and vizplotgui_tool.py commands, replace "tool"
with "gl" or "atomeye" or "pymol", depending on what visualization
package you have installed. We hope to add a VMD option soon.
with "gl" or "atomeye" or "pymol" or "vmd", depending on what
visualization package you have installed.
Note that for GL, you need to be able to run the Pizza.py GL tool,
which is included in the pizza sub-directory. See the "Pizza.py doc
@ -616,7 +616,7 @@ pages"_pizza for more info:
:link(pizza,http://www.sandia.gov/~sjplimp/pizza.html)
Note that for AtomEye, you need version 3, and their is a line in the
Note that for AtomEye, you need version 3, and there is a line in the
scripts that specifies the path and name of the executable. See the
AtomEye WWW pages "here"_atomeye or "here"_atomeye3 for more details:
@ -642,6 +642,10 @@ http://sourceforge.net/scm/?type=svn&group_id=4546 :pre
The latter link is to the open-source version.
Note that for VMD, you need a fairly current version (1.8.7 works for
me) and there are some lines in the pizza/vmd.py script for 4 PIZZA
variables that have to match the VMD installation on your system.
:line
See the python/README file for instructions on how to run them and the
@ -653,4 +657,4 @@ different visualization package options. Click to see larger images:
:image(JPG/screenshot_gl_small.jpg,JPG/screenshot_gl.jpg)
:image(JPG/screenshot_atomeye_small.jpg,JPG/screenshot_atomeye.jpg)
:image(JPG/screenshot_pymol_small.jpg,JPG/screenshot_pymol.jpg)
:image(JPG/screenshot_vmd_small.jpg,JPG/screenshot_vmd.jpg)

View File

@ -379,8 +379,8 @@ and barostatting to the system's potential energy as part of
</P>
<P>These fixes compute a global scalar and a global vector of quantities,
which can be accessed by various <A HREF = "Section_howto.html#4_15">output
commands</A>. The scalar values calculated by
this fix are "extensive"; the vector values are "intensive".
commands</A>. The scalar value calculated by
these fixes is "extensive"; the vector values are "intensive".
</P>
<P>The scalar is the cumulative energy change due to the fix.
</P>

View File

@ -371,8 +371,8 @@ and barostatting to the system's potential energy as part of
These fixes compute a global scalar and a global vector of quantities,
which can be accessed by various "output
commands"_Section_howto.html#4_15. The scalar values calculated by
this fix are "extensive"; the vector values are "intensive".
commands"_Section_howto.html#4_15. The scalar value calculated by
these fixes is "extensive"; the vector values are "intensive".
The scalar is the cumulative energy change due to the fix.
@ -433,7 +433,6 @@ limiting factor for numerical stability. Both
factorizations are time-reversible and can be shown to preserve the phase
space measure of the underlying non-Hamiltonian equations of motion.
[Restrictions:]
Non-periodic dimensions cannot be barostatted. {Z}, {xz}, and {yz},

View File

@ -272,21 +272,28 @@ for info on how to re-specify a fix in an input script that reads a
restart file, so that the operation of the fix continues in an
uninterrupted fashion.
</P>
<P>None of the <A HREF = "fix_modify.html">fix_modify</A> options are relevant to these
fixes.
<P>The <A HREF = "fix_modify.html">fix_modify</A> <I>energy</I> option is supported by the
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>These fixes compute a global array of values which can be accessed by
various <A HREF = "Section_howto.html#4_15">output commands</A>. The number of rows
in the array is equal to the number of rigid bodies. The number of
columns is 15. Thus for each rigid body, 12 values are stored: the
xyz coords of the center of mass (COM), the xyz components of the COM
velocity, the xyz components of the force acting on the COM, the xyz
components of the torque acting on the COM, and the xyz image flags of
the COM, which have the same meaning as image flags for atom positions
(see the "dump" command). The force and torque values in the array
are not affected by the <I>force</I> and <I>torque</I> keywords in the fix rigid
command; they reflect values before any changes are made by those
keywords.
<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
cumulative energy change due to the thermostatting the fix performs.
</P>
<P>All of these fixes compute a global array of values which can be
accessed by various <A HREF = "Section_howto.html#4_15">output commands</A>. The
number of rows in the array is equal to the number of rigid bodies.
The number of columns is 15. Thus for each rigid body, 12 values are
stored: the xyz coords of the center of mass (COM), the xyz components
of the COM velocity, the xyz components of the force acting on the
COM, the xyz components of the torque acting on the COM, and the xyz
image flags of the COM, which have the same meaning as image flags for
atom positions (see the "dump" command). The force and torque values
in the array are not affected by the <I>force</I> and <I>torque</I> keywords in
the fix rigid command; they reflect values before any changes are made
by those keywords.
</P>
<P>The ordering of the rigid bodies (by row in the array) is as follows.
For the <I>single</I> keyword there is just one rigid body. For the

View File

@ -261,21 +261,28 @@ for info on how to re-specify a fix in an input script that reads a
restart file, so that the operation of the fix continues in an
uninterrupted fashion.
None of the "fix_modify"_fix_modify.html options are relevant to these
fixes.
The "fix_modify"_fix_modify.html {energy} option is supported by the
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.
These fixes compute a global array of values which can be accessed by
various "output commands"_Section_howto.html#4_15. The number of rows
in the array is equal to the number of rigid bodies. The number of
columns is 15. Thus for each rigid body, 12 values are stored: the
xyz coords of the center of mass (COM), the xyz components of the COM
velocity, the xyz components of the force acting on the COM, the xyz
components of the torque acting on the COM, and the xyz image flags of
the COM, which have the same meaning as image flags for atom positions
(see the "dump" command). The force and torque values in the array
are not affected by the {force} and {torque} keywords in the fix rigid
command; they reflect values before any changes are made by those
keywords.
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
cumulative energy change due to the thermostatting the fix performs.
All of these fixes compute a global array of values which can be
accessed by various "output commands"_Section_howto.html#4_15. The
number of rows in the array is equal to the number of rigid bodies.
The number of columns is 15. Thus for each rigid body, 12 values are
stored: the xyz coords of the center of mass (COM), the xyz components
of the COM velocity, the xyz components of the force acting on the
COM, the xyz components of the torque acting on the COM, and the xyz
image flags of the COM, which have the same meaning as image flags for
atom positions (see the "dump" command). The force and torque values
in the array are not affected by the {force} and {torque} keywords in
the fix rigid command; they reflect values before any changes are made
by those keywords.
The ordering of the rigid bodies (by row in the array) is as follows.
For the {single} keyword there is just one rigid body. For the

View File

@ -52,11 +52,13 @@
group functions = count(group), mass(group), charge(group),
xcm(group,dim), vcm(group,dim), fcm(group,dim),
bound(group,xmin), gyration(group), ke(group),
angmom(group,dim),inertia(group,dimdim),omega(group,dim)
angmom(group,dim), torque(group,dim),
inertia(group,dimdim), omega(group,dim)
region functions = count(group,region), mass(group,region), charge(group,region),
xcm(group,dim,region), vcm(group,dim,region), fcm(group,dim,region),
bound(group,xmin,region), gyration(group,region), ke(group,reigon),
angmom(group,dim,region), inertia(group,dimdim,region),omega(group,dim,region)
angmom(group,dim,region), torque(group,dim,region),
inertia(group,dimdim,region), omega(group,dim,region)
special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)
atom value = mass[i], type[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i]
atom vector = mass, type, x, y, z, vx, vy, vz, fx, fy, fz
@ -274,8 +276,8 @@ references to other variables.
<TR><TD >Thermo keywords</TD><TD > vol, pe, ebond, etc</TD></TR>
<TR><TD >Math operators</TD><TD > (), -x, x+y, x-y, x*y, x/y, x^y, x==y, x!=y, x<y, x<=y, x>y, x>=y, x&&y, x||y, !x</TD></TR>
<TR><TD >Math functions</TD><TD > sqrt(x), exp(x), ln(x), log(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)</TD></TR>
<TR><TD >Group functions</TD><TD > count(ID), mass(ID), charge(ID), xcm(ID,dim), vcm(ID,dim), fcm(ID,dim), bound(ID,dir), gyration(ID), ke(ID), angmom(ID,dim), inertia(ID,dimdim), omega(ID,dim)</TD></TR>
<TR><TD >Region functions</TD><TD > count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR)</TD></TR>
<TR><TD >Group functions</TD><TD > count(ID), mass(ID), charge(ID), xcm(ID,dim), vcm(ID,dim), fcm(ID,dim), bound(ID,dir), gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), inertia(ID,dimdim), omega(ID,dim)</TD></TR>
<TR><TD >Region functions</TD><TD > count(ID,IDR), mass(ID,IDR), charge(ID,IDR), xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), angmom(ID,dim,IDR), torque(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR)</TD></TR>
<TR><TD >Special functions</TD><TD > sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)</TD></TR>
<TR><TD >Atom values</TD><TD > mass[i], type[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i]</TD></TR>
<TR><TD >Atom vectors</TD><TD > mass, type, x, y, z, vx, vy, vz, fx, fy, fz</TD></TR>
@ -476,10 +478,12 @@ the min/max of a particular coordinate for all atoms in the group.
Gyration() computes the radius-of-gyration of the group of atoms. See
the <A HREF = "compute_gyration.html">compute gyration</A> command for a definition
of the formula. Angmom() returns components of the angular momentum
of the group of atoms around its center of mass. Inertia() returns
one of 6 components of the inertia tensor of the group of atoms around
its center of mass. Omega() returns components of the angular
velocity of the group of atoms around its center of mass.
of the group of atoms around its center of mass. Torque() returns
components of the torque on the group of atoms around its center of
mass, based on current forces on the atoms. Inertia() returns one of
6 components of the inertia tensor of the group of atoms around its
center of mass. Omega() returns components of the angular velocity of
the group of atoms around its center of mass.
</P>
<P>Region functions are specified exactly the same way as group functions
except they take an extra argument which is the region ID. The

View File

@ -47,11 +47,13 @@ style = {delete} or {index} or {loop} or {world} or {universe} or {uloop} or {st
group functions = count(group), mass(group), charge(group),
xcm(group,dim), vcm(group,dim), fcm(group,dim),
bound(group,xmin), gyration(group), ke(group),
angmom(group,dim),inertia(group,dimdim),omega(group,dim)
angmom(group,dim), torque(group,dim),
inertia(group,dimdim), omega(group,dim)
region functions = count(group,region), mass(group,region), charge(group,region),
xcm(group,dim,region), vcm(group,dim,region), fcm(group,dim,region),
bound(group,xmin,region), gyration(group,region), ke(group,reigon),
angmom(group,dim,region), inertia(group,dimdim,region),omega(group,dim,region)
angmom(group,dim,region), torque(group,dim,region),
inertia(group,dimdim,region), omega(group,dim,region)
special functions = sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)
atom value = mass\[i\], type\[i\], x\[i\], y\[i\], z\[i\], vx\[i\], vy\[i\], vz\[i\], fx\[i\], fy\[i\], fz\[i\]
atom vector = mass, type, x, y, z, vx, vy, vz, fx, fy, fz
@ -269,12 +271,13 @@ Math operators: (), -x, x+y, x-y, x*y, x/y, x^y, x==y, x!=y, x<y, x<=y, x>y, x>=
Math functions: sqrt(x), exp(x), ln(x), log(x), sin(x), cos(x), tan(x), asin(x), acos(x), atan(x), atan2(y,x), random(x,y,z), normal(x,y,z), ceil(x), floor(x), round(x), ramp(x,y), stagger(x,y), logfreq(x,y,z), vdisplace(x,y), swiggle(x,y,z), cwiggle(x,y,z)
Group functions: count(ID), mass(ID), charge(ID), xcm(ID,dim), \
vcm(ID,dim), fcm(ID,dim), bound(ID,dir), \
gyration(ID), ke(ID), angmom(ID,dim), \
gyration(ID), ke(ID), angmom(ID,dim), torque(ID,dim), \
inertia(ID,dimdim), omega(ID,dim)
Region functions: count(ID,IDR), mass(ID,IDR), charge(ID,IDR), \
xcm(ID,dim,IDR), vcm(ID,dim,IDR), fcm(ID,dim,IDR), \
bound(ID,dir,IDR), gyration(ID,IDR), ke(ID,IDR), \
angmom(ID,dim,IDR), inertia(ID,dimdim,IDR), omega(ID,dim,IDR)
angmom(ID,dim,IDR), torque(ID,dim,IDR), \
inertia(ID,dimdim,IDR), omega(ID,dim,IDR)
Special functions: sum(x), min(x), max(x), ave(x), trap(x), gmask(x), rmask(x), grmask(x,y)
Atom values: mass\[i\], type\[i\], x\[i\], y\[i\], z\[i\], \
vx\[i\], vy\[i\], vz\[i\], fx\[i\], fy\[i\], fz\[i\]
@ -475,10 +478,12 @@ the min/max of a particular coordinate for all atoms in the group.
Gyration() computes the radius-of-gyration of the group of atoms. See
the "compute gyration"_compute_gyration.html command for a definition
of the formula. Angmom() returns components of the angular momentum
of the group of atoms around its center of mass. Inertia() returns
one of 6 components of the inertia tensor of the group of atoms around
its center of mass. Omega() returns components of the angular
velocity of the group of atoms around its center of mass.
of the group of atoms around its center of mass. Torque() returns
components of the torque on the group of atoms around its center of
mass, based on current forces on the atoms. Inertia() returns one of
6 components of the inertia tensor of the group of atoms around its
center of mass. Omega() returns components of the angular velocity of
the group of atoms around its center of mass.
Region functions are specified exactly the same way as group functions
except they take an extra argument which is the region ID. The

View File

@ -1,13 +1,20 @@
# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html
# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
#
# Minimalistic VMD embedding for Pizza.py
#
# This class will replace the VMD startup script,
# open a pipe to the executable, and feed it Tcl
# command lines one at a time.
#
# (c) 2010 Axel Kohlmeyer <akohlmey@gmail.com>
# Copyright (2005) Sandia Corporation. Under the terms of Contract
# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
# certain rights in this software. This software is distributed under
# the GNU General Public License.
oneline = "Control VMD from python."
# vmd tool
# Minimalistic VMD embedding for Pizza.py
# (c) 2010 Axel Kohlmeyer <akohlmey@gmail.com>
# This class will replace the VMD startup script,
# open a pipe to the executable,
# and feed it Tcl command lines one at a time
oneline = "Control VMD from python"
docstr = """
v = vmd() start up VMD
@ -30,6 +37,9 @@ v.enter() enter interactive shell
v.debug([True|False]) display generated VMD script commands?
"""
# History
# 11/10, Axel Kohlmeyer (Temple U): original version
# Imports and external programs
import types, os
@ -42,12 +52,14 @@ except: PIZZA_VMDDIR = "/usr/local/lib/vmd"
try: from DEFAULTS import PIZZA_VMDDEV
except: PIZZA_VMDDEV = "win"
try: from DEFAULTS import PIZZA_VMDARCH
except: PIZZA_VMDARCH = "LINUXAMD64"
except: PIZZA_VMDARCH = "LINUX"
try: import pexpect
except:
print "pexpect from http://pypi.python.org/pypi/pexpect/ is required for vmd module"
print "pexpect from http://pypi.python.org/pypi/pexpect", \
"is required for vmd tool"
raise
# Class definition
class vmd:

View File

@ -1,9 +1,9 @@
#!/usr/local/bin/python -i
# preceeding line should have path for Python on your machine
# viz.py
# Purpose: viz running LAMMPS simulation via GL tool in Pizza.py
# Syntax: viz.py in.lammps Nfreq Nsteps
# viz_atomeye.py
# Purpose: viz running LAMMPS simulation via AtomEye
# Syntax: viz_atomeye.py in.lammps Nfreq Nsteps
# in.lammps = LAMMPS input script
# Nfreq = dump and viz shapshot every this many steps
# Nsteps = run for this many steps
@ -11,14 +11,15 @@
import sys,os
# set this to point to AtomEye version 3 executable
ATOMEYE3 = "/home/sjplimp/tools/atomeye3/A3.i686-20060530"
# first line if want AtomEye output to screen, 2nd line to file
#ATOMEYE3 = "/home/sjplimp/tools/atomeye3/A3.i686-20060530"
ATOMEYE3 = "/home/sjplimp/tools/atomeye3/A3.i686-20060530 > atomeye.out"
# parse command line
argv = sys.argv
if len(argv) != 4:
print "Syntax: viz.py in.lammps Nfreq Nsteps"
print "Syntax: viz_atomeye.py in.lammps Nfreq Nsteps"
sys.exit()
infile = sys.argv[1]

View File

@ -1,9 +1,9 @@
#!/usr/local/bin/python -i
# preceeding line should have path for Python on your machine
# viz.py
# viz_gl.py
# Purpose: viz running LAMMPS simulation via GL tool in Pizza.py
# Syntax: viz.py in.lammps Nfreq Nsteps
# Syntax: viz_gl.py in.lammps Nfreq Nsteps
# in.lammps = LAMMPS input script
# Nfreq = dump and viz shapshot every this many steps
# Nsteps = run for this many steps
@ -15,7 +15,7 @@ sys.path.append("./pizza")
argv = sys.argv
if len(argv) != 4:
print "Syntax: viz.py in.lammps Nfreq Nsteps"
print "Syntax: viz_gl.py in.lammps Nfreq Nsteps"
sys.exit()
infile = sys.argv[1]

View File

@ -81,9 +81,9 @@ lmp.command("run 0 pre no post yes")
if me == 0:
v.flush()
# uncomment the following, if you want to work with the viz some more.
v('menu main on')
# print "type quit to terminate."
v.enter()
#v('menu main on')
#print "type quit to terminate."
#v.enter()
#v.stop()
# uncomment if running in parallel via Pypar

View File

@ -1,9 +1,9 @@
#!/usr/local/bin/python -i
# preceeding line should have path for Python on your machine
# vizplotgui.py
# Purpose: viz running LAMMPS simulation with plot and GUI
# Syntax: vizplotgui.py in.lammps Nfreq compute-ID
# vizplotgui_atomeye.py
# Purpose: viz running LAMMPS simulation via AtomEye with plot and GUI
# Syntax: vizplotgui_atomeye.py in.lammps Nfreq compute-ID
# in.lammps = LAMMPS input script
# Nfreq = plot data point and viz shapshot every this many steps
# compute-ID = ID of compute that calculates temperature
@ -15,6 +15,7 @@
import sys,os,time
sys.path.append("./pizza")
# set this to point to AtomEye version 3 executable
# first line if want AtomEye output to screen, 2nd line to file
#ATOMEYE3 = "/home/sjplimp/tools/atomeye3/A3.i686-20060530"
ATOMEYE3 = "/home/sjplimp/tools/atomeye3/A3.i686-20060530 > atomeye.out"
@ -49,7 +50,7 @@ def update(ntimestep):
argv = sys.argv
if len(argv) != 4:
print "Syntax: vizplotgui.py in.lammps Nfreq compute-ID"
print "Syntax: vizplotgui_atomeye.py in.lammps Nfreq compute-ID"
sys.exit()
infile = sys.argv[1]
@ -85,7 +86,7 @@ breakflag = 0
runflag = 0
temptarget = 1.0
# wrapper on GL window via Pizza.py gl tool
# wrapper on AtomEye
# just proc 0 handles reading of dump file and viz
if me == 0:

View File

@ -1,9 +1,9 @@
#!/usr/local/bin/python -i
# preceeding line should have path for Python on your machine
# vizplotgui.py
# Purpose: viz running LAMMPS simulation with plot and GUI
# Syntax: vizplotgui.py in.lammps Nfreq compute-ID
# vizplotgui_gl.py
# Purpose: viz running LAMMPS simulation via GL tool with plot and GUI
# Syntax: vizplotgui_gl.py in.lammps Nfreq compute-ID
# in.lammps = LAMMPS input script
# Nfreq = plot data point and viz shapshot every this many steps
# compute-ID = ID of compute that calculates temperature
@ -46,7 +46,7 @@ def update(ntimestep):
argv = sys.argv
if len(argv) != 4:
print "Syntax: vizplotgui.py in.lammps Nfreq compute-ID"
print "Syntax: vizplotgui_gl.py in.lammps Nfreq compute-ID"
sys.exit()
infile = sys.argv[1]

View File

@ -2,7 +2,7 @@
# preceeding line should have path for Python on your machine
# vizplotgui_pymol.py
# Purpose: viz running LAMMPS simulation with plot and GUI
# Purpose: viz running LAMMPS simulation via PyMol with plot and GUI
# Syntax: vizplotgui_pymol.py in.lammps Nfreq compute-ID
# in.lammps = LAMMPS input script
# Nfreq = plot data point and viz shapshot every this many steps
@ -48,7 +48,7 @@ def update(ntimestep):
argv = sys.argv
if len(argv) != 4:
print "Syntax: vizplotgui.py in.lammps Nfreq compute-ID"
print "Syntax: vizplotgui_pymol.py in.lammps Nfreq compute-ID"
sys.exit()
infile = sys.argv[1]

173
python/examples/vizplotgui_vmd.py Executable file
View File

@ -0,0 +1,173 @@
#!/usr/local/bin/python -i
# preceeding line should have path for Python on your machine
# vizplotgui_vmd.py
# Purpose: viz running LAMMPS simulation via VMD with plot and GUI
# Syntax: vizplotgui_vmd.py in.lammps Nfreq compute-ID
# in.lammps = LAMMPS input script
# Nfreq = plot data point and viz shapshot every this many steps
# compute-ID = ID of compute that calculates temperature
# (or any other scalar quantity)
# IMPORTANT: this script cannot yet be run in parallel via Pypar,
# because I can't seem to do a MPI-style broadcast in Pypar
import sys,time
sys.path.append("./pizza")
# methods called by GUI
def run():
global runflag
runflag = 1
def stop():
global runflag
runflag = 0
def settemp(value):
global temptarget
temptarget = slider.get()
def quit():
global breakflag
breakflag = 1
# method called by timestep loop every Nfreq steps
# read dump snapshot and viz it, update plot with compute value
def update(ntimestep):
d.next()
d.unscale()
p.single(ntimestep)
v.append('tmp.pdb','pdb')
value = lmp.extract_compute(compute,0,0)
xaxis.append(ntimestep)
yaxis.append(value)
gn.plot(xaxis,yaxis)
# parse command line
argv = sys.argv
if len(argv) != 4:
print "Syntax: vizplotgui_vmd.py in.lammps Nfreq compute-ID"
sys.exit()
infile = sys.argv[1]
nfreq = int(sys.argv[2])
compute = sys.argv[3]
me = 0
# uncomment if running in parallel via Pypar
#import pypar
#me = pypar.rank()
#nprocs = pypar.size()
from lammps import lammps
lmp = lammps()
# run infile all at once
# assumed to have no run command in it
# dump a file in native LAMMPS dump format for Pizza.py dump tool
lmp.file(infile)
lmp.command("thermo %d" % nfreq)
lmp.command("dump python all atom %d tmp.dump" % nfreq)
# initial 0-step run to generate initial 1-point plot, dump file, and image
lmp.command("run 0 pre yes post no")
value = lmp.extract_compute(compute,0,0)
ntimestep = 0
xaxis = [ntimestep]
yaxis = [value]
breakflag = 0
runflag = 0
temptarget = 1.0
# wrapper on VMD window via Pizza.py vmd tool
# just proc 0 handles reading of dump file and viz
if me == 0:
from vmd import vmd
v = vmd()
v('menu main off')
v.rep('VDW')
from dump import dump
from pdbfile import pdbfile
d = dump('tmp.dump',0)
p = pdbfile(d)
d.next()
d.unscale()
p.single(ntimestep)
v.new('tmp.pdb','pdb')
# display GUI with run/stop buttons and slider for temperature
if me == 0:
from Tkinter import *
tkroot = Tk()
tkroot.withdraw()
root = Toplevel(tkroot)
root.title("LAMMPS GUI")
frame = Frame(root)
Button(frame,text="Run",command=run).pack(side=LEFT)
Button(frame,text="Stop",command=stop).pack(side=LEFT)
slider = Scale(frame,from_=0.0,to=5.0,resolution=0.1,
orient=HORIZONTAL,label="Temperature")
slider.bind('<ButtonRelease-1>',settemp)
slider.set(temptarget)
slider.pack(side=LEFT)
Button(frame,text="Quit",command=quit).pack(side=RIGHT)
frame.pack()
tkroot.update()
# wrapper on GnuPlot via Pizza.py gnu tool
if me == 0:
from gnu import gnu
gn = gnu()
gn.plot(xaxis,yaxis)
gn.title(compute,"Timestep","Temperature")
# endless loop, checking status of GUI settings every Nfreq steps
# run with pre yes/no and post yes/no depending on go/stop status
# re-invoke fix langevin with new seed when temperature slider changes
# after re-invoke of fix langevin, run with pre yes
running = 0
temp = temptarget
seed = 12345
lmp.command("fix 2 all langevin %g %g 0.1 %d" % (temp,temp,seed))
while 1:
if me == 0: tkroot.update()
if temp != temptarget:
temp = temptarget
seed += me+1
lmp.command("fix 2 all langevin %g %g 0.1 12345" % (temp,temp))
running = 0
if runflag and running:
lmp.command("run %d pre no post no" % nfreq)
ntimestep += nfreq
if me == 0: update(ntimestep)
elif runflag and not running:
lmp.command("run %d pre yes post no" % nfreq)
ntimestep += nfreq
if me == 0: update(ntimestep)
elif not runflag and running:
lmp.command("run %d pre no post yes" % nfreq)
ntimestep += nfreq
if me == 0: update(ntimestep)
if breakflag: break
if runflag: running = 1
else: running = 0
time.sleep(0.01)
lmp.command("run 0 pre no post yes")
# uncomment if running in parallel via Pypar
#print "Proc %d out of %d procs has" % (me,nprocs), lmp
#pypar.finalize()

View File

@ -275,6 +275,8 @@ void PairCoulLong::init_style()
double PairCoulLong::init_one(int i, int j)
{
scale[j][i] = scale[i][j];
return cut_coul;
}

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (U Penn), akohlmey at cmm.chem.upenn.edu
Contributing author: Axel Kohlmeyer (Temple U), akohlmey at gmail.com
------------------------------------------------------------------------- */
#include "math.h"

View File

@ -132,6 +132,7 @@ void FixNEB::min_setup(int vflag)
void FixNEB::min_post_force(int vflag)
{
MPI_Status status;
MPI_Request request;
double vprev,vnext,vmax,vmin;
double delx,dely,delz;
double delta1[3],delta2[3];
@ -141,10 +142,21 @@ void FixNEB::min_post_force(int vflag)
veng = pe->compute_scalar();
if (ireplica < nreplica-1)
if (ireplica == 0)
MPI_Send(&veng,1,MPI_DOUBLE,procnext,0,uworld);
else if (ireplica == nreplica-1) {
MPI_Irecv(&vprev,1,MPI_DOUBLE,procprev,0,uworld,&request);
MPI_Wait(&request,&status);
} else
MPI_Sendrecv(&veng,1,MPI_DOUBLE,procnext,0,
&vprev,1,MPI_DOUBLE,procprev,0,uworld,&status);
if (ireplica > 0)
if (ireplica == 0) {
MPI_Irecv(&vnext,1,MPI_DOUBLE,procnext,0,uworld,&request);
MPI_Wait(&request,&status);
} else if (ireplica == nreplica-1)
MPI_Send(&veng,1,MPI_DOUBLE,procprev,0,uworld);
else
MPI_Sendrecv(&veng,1,MPI_DOUBLE,procprev,0,
&vnext,1,MPI_DOUBLE,procnext,0,uworld,&status);
@ -157,10 +169,21 @@ void FixNEB::min_post_force(int vflag)
int nlocal = atom->nlocal;
if (nlocal != natoms) error->one("Atom count changed in fix neb");
if (ireplica < nreplica-1)
if (ireplica == 0)
MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld);
else if (ireplica == nreplica-1) {
MPI_Irecv(xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&request);
MPI_Wait(&request,&status);
} else
MPI_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procnext,0,
xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&status);
if (ireplica > 0)
if (ireplica == 0) {
MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request);
MPI_Wait(&request,&status);
} else if (ireplica == nreplica-1)
MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld);
else
MPI_Sendrecv(x[0],3*nlocal,MPI_DOUBLE,procprev,0,
xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&status);

View File

@ -383,12 +383,13 @@ void NEB::print_status()
for (int i = 1; i < nreplica; i++)
rdist[i] = rdist[i-1] + all[i][1];
double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2];
if (endpt > 0.0)
for (int i = 1; i < nreplica; i++)
rdist[i] /= endpt;
// look up GradV for the initial, final, and climbing replicas
// These are identical to fnorm2, but better to be safe we
// take them straight from fix_neb
// these should be identical to fnorm2
// but to be safe take them straight from fix neb
double gradvnorm0, gradvnorm1, gradvnormc;

View File

@ -13,7 +13,7 @@
/* ----------------------------------------------------------------------
Common parameters for the CMM coarse grained MD potentials.
Contributing author: Axel Kohlmeyer <akohlmey@gmail.com>
Contributing author: Axel Kohlmeyer (Temple U) <akohlmey@gmail.com>
------------------------------------------------------------------------- */
#ifndef LMP_CG_CMM_PARMS_H

View File

@ -1416,6 +1416,13 @@ double FixNH::compute_vector(int n)
/* ---------------------------------------------------------------------- */
void FixNH::reset_target(double t_new)
{
t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
void FixNH::reset_dt()
{
dtv = update->dt;

View File

@ -34,6 +34,7 @@ class FixNH : public Fix {
void write_restart(FILE *);
void restart(char *);
int modify_param(int, char **);
void reset_target(double);
void reset_dt();
protected:

View File

@ -28,7 +28,7 @@ class FixRigid : public Fix {
public:
FixRigid(class LAMMPS *, int, char **);
virtual ~FixRigid();
int setmask();
virtual int setmask();
virtual void init();
virtual void setup(int);
virtual void initial_integrate(int);

View File

@ -42,7 +42,9 @@ FixRigidNVT::FixRigidNVT(LAMMPS *lmp, int narg, char **arg) :
{
// other settings are made by FixRigid parent
scalar_flag = 1;
restart_global = 1;
extscalar = 1;
// error checking
// convert input period to frequency
@ -86,6 +88,20 @@ FixRigidNVT::~FixRigidNVT()
/* ---------------------------------------------------------------------- */
int FixRigidNVT::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= PRE_NEIGHBOR;
mask |= THERMO_ENERGY;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixRigidNVT::init()
{
FixRigid::init();
@ -105,6 +121,8 @@ void FixRigidNVT::init()
for (int k = 0; k < domain->dimension; k++)
if (fabs(inertia[ibody][k]) < 1e-6) nf_r--;
// see Table 1 in Kamberaj et al
if (t_order == 3) {
w[0] = 1.0 / (2.0 - pow(2.0, 1.0/3.0));
w[1] = 1.0 - 2.0*w[0];
@ -200,12 +218,12 @@ void FixRigidNVT::initial_integrate(int vflag)
xcm[ibody][1] += dtv * vcm[ibody][1];
xcm[ibody][2] += dtv * vcm[ibody][2];
// step 1.3 - apply torque (body coords) to quaternion momentum
torque[ibody][0] *= tflag[ibody][0];
torque[ibody][1] *= tflag[ibody][1];
torque[ibody][2] *= tflag[ibody][2];
// step 1.3 - apply torque (body coords) to quaternion momentum
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody);
quatvec(quat[ibody],tbody,fquat);
@ -219,7 +237,7 @@ void FixRigidNVT::initial_integrate(int vflag)
conjqm[ibody][2] *= scale_r;
conjqm[ibody][3] *= scale_r;
// step 1.4 to 1.13 - use no_squish rotate to update p and q
// step 1.4 to 1.8 - use no_squish rotate to update p (i.e. conjqm) and q
no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
@ -227,7 +245,7 @@ void FixRigidNVT::initial_integrate(int vflag)
no_squish_rotate(2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
// update the exyz_space
// update the exyz_space from new quaternion
// transform p back to angmom
// update angular velocity
@ -271,7 +289,7 @@ void FixRigidNVT::final_integrate()
double tmp,scale_t,scale_r,akin_t,akin_r;
double dtfm,xy,xz,yz;
// intialize velocity scale for translation and rotation
// compute velocity scales for translation and rotation
tmp = -1.0 * dtq * eta_dot_t[0];
scale_t = exp(tmp);
@ -360,7 +378,7 @@ void FixRigidNVT::final_integrate()
torque[ibody][1] = all[ibody][4];
torque[ibody][2] = all[ibody][5];
// update vcm by 1/2 step
// 2.5-2.6 update vcm by 1/2 step
dtfm = dtf / masstotal[ibody];
vcm[ibody][0] *= scale_t;
@ -370,22 +388,31 @@ void FixRigidNVT::final_integrate()
vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
// update conjqm, then transform to angmom, set velocity again
// virial is already setup from initial_integrate
// 2.1-2.4 update conjqm, angular momentum and angular velocity
// apply body torque flags
torque[ibody][0] *= tflag[ibody][0];
torque[ibody][1] *= tflag[ibody][1];
torque[ibody][2] *= tflag[ibody][2];
// convert torque to the body frame
matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody],
torque[ibody],tbody);
// compute "force" for quaternion
quatvec(quat[ibody],tbody,fquat);
// update the conjugate quaternion momentum (conjqm)
conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0];
conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1];
conjqm[ibody][2] = scale_r * conjqm[ibody][2] + dtf2 * fquat[2];
conjqm[ibody][3] = scale_r * conjqm[ibody][3] + dtf2 * fquat[3];
// compute angular momentum in the body frame then convert to the space-fixed frame
invquatvec(quat[ibody],conjqm[ibody],mbody);
matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody],
mbody,angmom[ibody]);
@ -394,6 +421,8 @@ void FixRigidNVT::final_integrate()
angmom[ibody][1] *= 0.5;
angmom[ibody][2] *= 0.5;
// compute new angular velocity
omega_from_angmom(angmom[ibody],ex_space[ibody],ey_space[ibody],
ez_space[ibody],inertia[ibody],omega[ibody]);
}
@ -535,6 +564,76 @@ void FixRigidNVT::write_restart(FILE *fp)
delete list;
}
/* ----------------------------------------------------------------------
compute kinetic energy in the extended Hamiltonian
conserved quantity = sum of returned energy and potential energy
-----------------------------------------------------------------------*/
double FixRigidNVT::compute_scalar()
{
int i,k,ibody;
double kt = boltz * t_target;
double energy,ke_t,ke_q,tmp,Pkq[4];
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
// translational kinetic energy
ke_t = 0.0;
for (ibody = 0; ibody < nbody; ibody++)
ke_t += 0.5 * masstotal[ibody] * (vcm[ibody][0]*vcm[ibody][0] +
vcm[ibody][1]*vcm[ibody][1] +
vcm[ibody][2]*vcm[ibody][2]);
// rotational kinetic energy
ke_q = 0.0;
for (ibody = 0; ibody < nbody; ibody++) {
for (k = 1; k < 4; k++) {
if (k == 1) {
Pkq[0] = -quat[ibody][1];
Pkq[1] = quat[ibody][0];
Pkq[2] = quat[ibody][3];
Pkq[3] = -quat[ibody][2];
} else if (k == 2) {
Pkq[0] = -quat[ibody][2];
Pkq[1] = -quat[ibody][3];
Pkq[2] = quat[ibody][0];
Pkq[3] = quat[ibody][1];
} else if (k == 3) {
Pkq[0] = -quat[ibody][3];
Pkq[1] = quat[ibody][2];
Pkq[2] = -quat[ibody][1];
Pkq[3] = quat[ibody][0];
}
tmp = conjqm[ibody][0]*Pkq[0] + conjqm[ibody][1]*Pkq[1] +
conjqm[ibody][2]*Pkq[2] + conjqm[ibody][3]*Pkq[3];
tmp *= tmp;
if (fabs(inertia[ibody][k-1]) < 1e-6) tmp = 0.0;
else tmp /= (8.0 * inertia[ibody][k-1]);
ke_q += tmp;
}
}
energy = ke_t + ke_q;
// thermostat chain energy: from equation 12 in Kameraj et al (JCP 2005)
energy += kt * (nf_t * eta_t[0] + nf_r * eta_r[0]);
for (i = 1; i < t_chain; i++)
energy += kt * (eta_t[i] + eta_r[i]);
for (i = 0; i < t_chain; i++) {
energy += 0.5 * q_t[i] * (eta_dot_t[i] * eta_dot_t[i]);
energy += 0.5 * q_r[i] * (eta_dot_r[i] * eta_dot_r[i]);
}
return energy;
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */

View File

@ -28,10 +28,12 @@ class FixRigidNVT : public FixRigid {
public:
FixRigidNVT(class LAMMPS *, int, char **);
virtual ~FixRigidNVT();
virtual int setmask();
virtual void init();
virtual void setup(int);
virtual void initial_integrate(int);
virtual void final_integrate();
virtual double compute_scalar();
virtual void write_restart(FILE *);
virtual void restart(char *);
virtual void reset_target(double);

View File

@ -1240,6 +1240,87 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
MPI_Allreduce(p,lmom,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute the torque T (tq) on group
around center-of-mass cm
must unwrap atoms to compute T correctly
------------------------------------------------------------------------- */
void Group::torque(int igroup, double *cm, double *tq)
{
int groupbit = bitmask[igroup];
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0];
}
MPI_Allreduce(tlocal,tq,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute the torque T (tq) on group of atoms in region
around center-of-mass cm
must unwrap atoms to compute T correctly
------------------------------------------------------------------------- */
void Group::torque(int igroup, double *cm, double *tq, int iregion)
{
int groupbit = bitmask[igroup];
Region *region = domain->regions[iregion];
double **x = atom->x;
double **f = atom->f;
int *mask = atom->mask;
int *image = atom->image;
int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0];
}
MPI_Allreduce(tlocal,tq,3,MPI_DOUBLE,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
compute moment of inertia tensor around center-of-mass cm of group
must unwrap atoms to compute itensor correctly

View File

@ -54,6 +54,8 @@ class Group : protected Pointers {
double gyration(int, double, double *, int);
void angmom(int, double *, double *); // angular momentum of group
void angmom(int, double *, double *, int);
void torque(int, double *, double *); // torque on group
void torque(int, double *, double *, int);
void inertia(int, double *, double [3][3]); // inertia tensor
void inertia(int, double *, double [3][3], int);
void omega(double *, double [3][3], double *); // angular velocity

View File

@ -219,6 +219,8 @@ double PairCoulCut::init_one(int i, int j)
if (setflag[i][j] == 0)
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
scale[j][i] = scale[i][j];
return cut[i][j];
}

View File

@ -2423,8 +2423,8 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
customize by adding a group function with optional region arg:
count(group),mass(group),charge(group),
xcm(group,dim),vcm(group,dim),fcm(group,dim),
bound(group,xmin),gyration(group),ke(group),angmom(group),
inertia(group,dim),omega(group,dim)
bound(group,xmin),gyration(group),ke(group),angmom(group,dim),
torque(group,dim),inertia(group,dim),omega(group,dim)
------------------------------------------------------------------------- */
int Variable::group_function(char *word, char *contents, Tree **tree,
@ -2438,7 +2438,8 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
strcmp(word,"vcm") && strcmp(word,"fcm") &&
strcmp(word,"bound") && strcmp(word,"gyration") &&
strcmp(word,"ke") && strcmp(word,"angmom") &&
strcmp(word,"inertia") && strcmp(word,"omega"))
strcmp(word,"torque") && strcmp(word,"inertia") &&
strcmp(word,"omega"))
return 0;
// parse contents for arg1,arg2,arg3 separated by commas
@ -2588,6 +2589,24 @@ int Variable::group_function(char *word, char *contents, Tree **tree,
else if (strcmp(arg2,"z") == 0) value = lmom[2];
else error->all("Invalid group function in variable formula");
} else if (strcmp(word,"torque") == 0) {
atom->check_mass();
double xcm[3],tq[3];
if (narg == 2) {
double masstotal = group->mass(igroup);
group->xcm(igroup,masstotal,xcm);
group->torque(igroup,xcm,tq);
} else if (narg == 3) {
int iregion = region_function(arg3);
double masstotal = group->mass(igroup,iregion);
group->xcm(igroup,masstotal,xcm,iregion);
group->torque(igroup,xcm,tq,iregion);
} else error->all("Invalid group function in variable formula");
if (strcmp(arg2,"x") == 0) value = tq[0];
else if (strcmp(arg2,"y") == 0) value = tq[1];
else if (strcmp(arg2,"z") == 0) value = tq[2];
else error->all("Invalid group function in variable formula");
} else if (strcmp(word,"inertia") == 0) {
atom->check_mass();
double xcm[3],inertia[3][3];

View File

@ -1 +1 @@
#define LAMMPS_VERSION "23 Nov 2010"
#define LAMMPS_VERSION "1 Dec 2010"

View File

@ -75,7 +75,7 @@ class Data {
char *atom_style;
int style_angle,style_atomic,style_bond,style_charge,style_dipole;
int style_dpd,style_ellipsoid,style_full,style_granular;
int style_ellipsoid,style_full,style_granular;
int style_hybrid,style_molecular,style_peri;
int natoms,nbonds,nangles,ndihedrals,nimpropers;
@ -204,7 +204,6 @@ class Data {
void write_atom_bond(FILE *, int, int, int, int);
void write_atom_charge(FILE *, int, int, int, int);
void write_atom_dipole(FILE *, int, int, int, int);
void write_atom_dpd(FILE *, int, int, int, int);
void write_atom_ellipsoid(FILE *, int, int, int, int);
void write_atom_full(FILE *, int, int, int, int);
void write_atom_granular(FILE *, int, int, int, int);
@ -216,7 +215,6 @@ class Data {
void write_atom_bond_extra(FILE *, int);
void write_atom_charge_extra(FILE *, int);
void write_atom_dipole_extra(FILE *, int);
void write_atom_dpd_extra(FILE *, int);
void write_atom_ellipsoid_extra(FILE *, int);
void write_atom_full_extra(FILE *, int);
void write_atom_granular_extra(FILE *, int);
@ -228,7 +226,6 @@ class Data {
void write_vel_bond(FILE *, int);
void write_vel_charge(FILE *, int);
void write_vel_dipole(FILE *, int);
void write_vel_dpd(FILE *, int);
void write_vel_ellipsoid(FILE *, int);
void write_vel_full(FILE *, int);
void write_vel_granular(FILE *, int);
@ -240,7 +237,6 @@ class Data {
void write_vel_bond_extra(FILE *, int);
void write_vel_charge_extra(FILE *, int);
void write_vel_dipole_extra(FILE *, int);
void write_vel_dpd_extra(FILE *, int);
void write_vel_ellipsoid_extra(FILE *, int);
void write_vel_full_extra(FILE *, int);
void write_vel_granular_extra(FILE *, int);
@ -270,7 +266,6 @@ void allocate_atomic(Data &data);
void allocate_bond(Data &data);
void allocate_charge(Data &data);
void allocate_dipole(Data &data);
void allocate_dpd(Data &data);
void allocate_ellipsoid(Data &data);
void allocate_full(Data &data);
void allocate_granular(Data &data);
@ -282,7 +277,6 @@ int atom_atomic(double *, Data &, int);
int atom_bond(double *, Data &, int);
int atom_charge(double *, Data &, int);
int atom_dipole(double *, Data &, int);
int atom_dpd(double *, Data &, int);
int atom_ellipsoid(double *, Data &, int);
int atom_full(double *, Data &, int);
int atom_granular(double *, Data &, int);
@ -464,7 +458,7 @@ void header(FILE *fp, Data &data)
else if (flag == ATOM_STYLE) {
data.style_angle = data.style_atomic = data.style_bond =
data.style_charge = data.style_dipole = data.style_dpd =
data.style_charge = data.style_dipole =
data.style_ellipsoid = data.style_full = data.style_granular =
data.style_hybrid = data.style_molecular = data.style_peri = 0;
@ -538,7 +532,6 @@ void set_style(char *name, Data &data, int flag)
else if (strcmp(name,"bond") == 0) data.style_bond = flag;
else if (strcmp(name,"charge") == 0) data.style_charge = flag;
else if (strcmp(name,"dipole") == 0) data.style_dipole = flag;
else if (strcmp(name,"dpd") == 0) data.style_dpd = flag;
else if (strcmp(name,"ellipsoid") == 0) data.style_ellipsoid = flag;
else if (strcmp(name,"full") == 0) data.style_full = flag;
else if (strcmp(name,"granular") == 0) data.style_granular = flag;
@ -712,7 +705,6 @@ int atom(double *buf, Data &data)
if (data.style_bond) allocate_bond(data);
if (data.style_charge) allocate_charge(data);
if (data.style_dipole) allocate_dipole(data);
if (data.style_dpd) allocate_dpd(data);
if (data.style_ellipsoid) allocate_ellipsoid(data);
if (data.style_full) allocate_full(data);
if (data.style_granular) allocate_granular(data);
@ -735,7 +727,6 @@ int atom(double *buf, Data &data)
if (k == data.style_bond) m += atom_bond(&buf[m],data,iatoms);
if (k == data.style_charge) m += atom_charge(&buf[m],data,iatoms);
if (k == data.style_dipole) m += atom_dipole(&buf[m],data,iatoms);
if (k == data.style_dpd) m += atom_dpd(&buf[m],data,iatoms);
if (k == data.style_ellipsoid) m += atom_ellipsoid(&buf[m],data,iatoms);
if (k == data.style_full) m += atom_full(&buf[m],data,iatoms);
if (k == data.style_granular) m += atom_granular(&buf[m],data,iatoms);
@ -892,23 +883,6 @@ int atom_dipole(double *buf, Data &data, int iatoms)
return m;
}
int atom_dpd(double *buf, Data &data, int iatoms)
{
int m = 1;
data.x[iatoms] = buf[m++];
data.y[iatoms] = buf[m++];
data.z[iatoms] = buf[m++];
data.tag[iatoms] = static_cast<int> (buf[m++]);
data.type[iatoms] = static_cast<int> (buf[m++]);
data.mask[iatoms] = static_cast<int> (buf[m++]);
data.image[iatoms] = static_cast<int> (buf[m++]);
data.vx[iatoms] = buf[m++];
data.vy[iatoms] = buf[m++];
data.vz[iatoms] = buf[m++];
return m;
}
int atom_ellipsoid(double *buf, Data &data, int iatoms)
{
int m = 1;
@ -1187,8 +1161,6 @@ void allocate_dipole(Data &data)
data.muz = new double[data.natoms];
}
void allocate_dpd(Data &data) {}
void allocate_full(Data &data)
{
data.q = new double[data.natoms];
@ -3362,7 +3334,6 @@ void Data::write(FILE *fp, FILE *fp2)
if (style_bond) write_atom_bond(fp,i,ix,iy,iz);
if (style_charge) write_atom_charge(fp,i,ix,iy,iz);
if (style_dipole) write_atom_dipole(fp,i,ix,iy,iz);
if (style_dpd) write_atom_dpd(fp,i,ix,iy,iz);
if (style_ellipsoid) write_atom_ellipsoid(fp,i,ix,iy,iz);
if (style_full) write_atom_full(fp,i,ix,iy,iz);
if (style_granular) write_atom_granular(fp,i,ix,iy,iz);
@ -3379,7 +3350,6 @@ void Data::write(FILE *fp, FILE *fp2)
if (k == style_bond) write_atom_bond_extra(fp,i);
if (k == style_charge) write_atom_charge_extra(fp,i);
if (k == style_dipole) write_atom_dipole_extra(fp,i);
if (k == style_dpd) write_atom_dpd_extra(fp,i);
if (k == style_ellipsoid) write_atom_ellipsoid_extra(fp,i);
if (k == style_full) write_atom_full_extra(fp,i);
if (k == style_granular) write_atom_granular_extra(fp,i);
@ -3401,7 +3371,6 @@ void Data::write(FILE *fp, FILE *fp2)
if (style_bond) write_vel_bond(fp,i);
if (style_charge) write_vel_charge(fp,i);
if (style_dipole) write_vel_dipole(fp,i);
if (style_dpd) write_vel_dpd(fp,i);
if (style_ellipsoid) write_vel_ellipsoid(fp,i);
if (style_full) write_vel_full(fp,i);
if (style_granular) write_vel_granular(fp,i);
@ -3417,7 +3386,6 @@ void Data::write(FILE *fp, FILE *fp2)
if (k == style_bond) write_vel_bond_extra(fp,i);
if (k == style_charge) write_vel_charge_extra(fp,i);
if (k == style_dipole) write_vel_dipole_extra(fp,i);
if (k == style_dpd) write_vel_dpd_extra(fp,i);
if (k == style_ellipsoid) write_vel_ellipsoid_extra(fp,i);
if (k == style_full) write_vel_full_extra(fp,i);
if (k == style_granular) write_vel_granular_extra(fp,i);
@ -3494,12 +3462,6 @@ void Data::write_atom_dipole(FILE *fp, int i, int ix, int iy, int iz)
tag[i],type[i],q[i],x[i],y[i],z[i],mux[i],muy[i],muz[i],ix,iy,iz);
}
void Data::write_atom_dpd(FILE *fp, int i, int ix, int iy, int iz)
{
fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %d %d %d",
tag[i],type[i],x[i],y[i],z[i],ix,iy,iz);
}
void Data::write_atom_ellipsoid(FILE *fp, int i, int ix, int iy, int iz)
{
fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
@ -3558,8 +3520,6 @@ void Data::write_atom_dipole_extra(FILE *fp, int i)
fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e",q[i],mux[i],muy[i],muz[i]);
}
void Data::write_atom_dpd_extra(FILE *fp, int i) {}
void Data::write_atom_ellipsoid_extra(FILE *fp, int i)
{
fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e",quatw[i],quati[i],quatj[i],quatk[i]);
@ -3615,11 +3575,6 @@ void Data::write_vel_dipole(FILE *fp, int i)
fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
}
void Data::write_vel_dpd(FILE *fp, int i)
{
fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
}
void Data::write_vel_ellipsoid(FILE *fp, int i)
{
fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
@ -3657,7 +3612,6 @@ void Data::write_vel_atomic_extra(FILE *fp, int i) {}
void Data::write_vel_bond_extra(FILE *fp, int i) {}
void Data::write_vel_charge_extra(FILE *fp, int i) {}
void Data::write_vel_dipole_extra(FILE *fp, int i) {}
void Data::write_vel_dpd_extra(FILE *fp, int i) {}
void Data::write_vel_ellipsoid_extra(FILE *fp, int i)
{

View File

@ -18,3 +18,14 @@ in the syntax file (lammps.vim). You can easily add new ones.
Gerolf Ziegenhain <gerolf@ziegenhain.com> 2007
---------------
updated by Sam Bateman, 11/2010
Sam Bateman
Naval Research Laboratory
Code 7434
1005 Balch Blvd.
Stennis Space Center, MS 39529
Phone: (228) 688-4328
Email: sam.bateman@nrlssc.navy.mil

View File

@ -1,20 +1,20 @@
" Vim syntax file
" Language: Lammps Simulation Script File
" Maintainer: Gerolf Ziegenhain <gerolf@ziegenhain.com>
" Updates: Axel Kohlmeyer <akohlmey@gmail.com>
" Latest Revision: 2010-06-06
" Updates: Axel Kohlmeyer <akohlmey@gmail.com>, Sam Bateman <sam.bateman@nrlssc.navy.mil>
" Latest Revision: 2010-11-25
syn clear
syn keyword lammpsOutput log write_restart dump undump thermo thermo_modify thermo_style print
syn keyword lammpsRead include read read_restart read_data
syn keyword lammpsLattice boundary units atom_style lattice region create_box create_atoms dielectric
syn keyword lammpsLattice delete_atoms change_box dimension newton replicate
syn keyword lammpsParticle pair_coeff pair_style pair_modify mass angle_coeff angle_style atom_modify
syn keyword lammpsParticle atom_style bond_coeff bond_style delete_bonds kspace_style kspace_modify
syn keyword lammpsParticle dihedral_style dihedral_coeff improper_style improper_coeff
syn keyword lammpsSetup min_style fix_modify run_style timestep neighbor fix unfix communicate
syn keyword lammpsSetup neigh_modify reset_timestep velocity
syn keyword lammpsLattice delete_atoms change_box dimension replicate
syn keyword lammpsParticle pair_coeff pair_style pair_modify mass velocity angle_coeff angle_style
syn keyword lammpsParticle atom_modify atom_style bond_coeff bond_style delete_bonds kspace_style
syn keyword lammpsParticle kspace_modify dihedral_style dihedral_coeff improper_style improper_coeff
syn keyword lammpsSetup min_style fix_modify run_style timestep neighbor neigh_modify fix unfix
syn keyword lammpsSetup communicate newton nthreads processors reset_timestep
syn keyword lammpsRun minimize run
syn keyword lammpsDefine variable group
@ -22,7 +22,9 @@ syn keyword lammpsRepeat jump next loop
syn keyword lammpsOperator equal add sub mult div
syn keyword lammpsConditional if then else
syn keyword lammpsConditional if then elif else
syn keyword lammpsSpecial EDGE NULL
syn region lammpsString start=+'+ end=+'+ oneline
syn region lammpsString start=+"+ end=+"+ oneline
@ -34,7 +36,7 @@ syn match lammpsFloat "\<[0-9]\+[edED][-+]\=[0-9]\+[ij]\=\>"
syn match lammpsComment "#.*$"
syn match lammpsVariable "\$\({[a-zA-Z0-9]\+}\)"
syn match lammpsVariable "\$\({[a-zA-Z0-9_]\+}\)"
syn match lammpsVariable "\$[A-Za-z]"
if !exists("did_lammps_syntax_inits")
@ -55,6 +57,7 @@ if !exists("did_lammps_syntax_inits")
hi link lammpsVariable Identifier
hi link lammpsConditional Conditional
hi link lammpsOperator Operator
hi link lammpsSpecial Number
endif
let b:current_syntax = "lammps"