diff --git a/bench-accel/bench_zbl/in.zbl b/bench-accel/bench_zbl/in.zbl
new file mode 100644
index 000000000..7fb27451d
--- /dev/null
+++ b/bench-accel/bench_zbl/in.zbl
@@ -0,0 +1,24 @@
+# 3d Lennard-Jones melt
+
+units metal
+atom_style atomic
+
+lattice fcc 4.0
+region box block 0 20 0 20 0 20
+create_box 1 box
+create_atoms 1 box
+mass 1 1.0
+
+velocity all create 300.0 87287 loop geom
+
+pair_style zbl 3.0 4.0
+pair_coeff 1 1 14.0
+
+neighbor 2.5 bin
+neigh_modify delay 5 every 1
+
+fix 1 all nve
+
+thermo 10
+
+run 100
diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt
index af357f094..c08fa47c1 100644
--- a/doc/Section_commands.txt
+++ b/doc/Section_commands.txt
@@ -1,1186 +1,1191 @@
<"Previous Section"_Section_start.html - "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc - "Next Section"_Section_packages.html :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
3. Commands :h3
This section describes how a LAMMPS input script is formatted and the
input script commands used to define a LAMMPS simulation.
3.1 "LAMMPS input script"_#cmd_1
3.2 "Parsing rules"_#cmd_2
3.3 "Input script structure"_#cmd_3
3.4 "Commands listed by category"_#cmd_4
3.5 "Commands listed alphabetically"_#cmd_5 :all(b)
:line
:line
3.1 LAMMPS input script :link(cmd_1),h4
LAMMPS executes by reading commands from a input script (text file),
one line at a time. When the input script ends, LAMMPS exits. Each
command causes LAMMPS to take some action. It may set an internal
variable, read in a file, or run a simulation. Most commands have
default settings, which means you only need to use the command if you
wish to change the default.
In many cases, the ordering of commands in an input script is not
important. However the following rules apply:
(1) LAMMPS does not read your entire input script and then perform a
simulation with all the settings. Rather, the input script is read
one line at a time and each command takes effect when it is read.
Thus this sequence of commands:
timestep 0.5
run 100
run 100 :pre
does something different than this sequence:
run 100
timestep 0.5
run 100 :pre
In the first case, the specified timestep (0.5 fmsec) is used for two
simulations of 100 timesteps each. In the 2nd case, the default
timestep (1.0 fmsec) is used for the 1st 100 step simulation and a 0.5
fmsec timestep is used for the 2nd one.
(2) Some commands are only valid when they follow other commands. For
example you cannot set the temperature of a group of atoms until atoms
have been defined and a group command is used to define which atoms
belong to the group.
(3) Sometimes command B will use values that can be set by command A.
This means command A must precede command B in the input script if it
is to have the desired effect. For example, the
"read_data"_read_data.html command initializes the system by setting
up the simulation box and assigning atoms to processors. If default
values are not desired, the "processors"_processors.html and
"boundary"_boundary.html commands need to be used before read_data to
tell LAMMPS how to map processors to the simulation box.
Many input script errors are detected by LAMMPS and an ERROR or
WARNING message is printed. "This section"_Section_errors.html gives
more information on what errors mean. The documentation for each
command lists restrictions on how the command can be used.
:line
3.2 Parsing rules :link(cmd_2),h4
Each non-blank line in the input script is treated as a command.
LAMMPS commands are case sensitive. Command names are lower-case, as
are specified command arguments. Upper case letters may be used in
file names or user-chosen ID strings.
Here is how each line in the input script is parsed by LAMMPS:
(1) If the last printable character on the line is a "&" character
(with no surrounding quotes), the command is assumed to continue on
the next line. The next line is concatenated to the previous line by
removing the "&" character and newline. This allows long commands to
be continued across two or more lines.
(2) All characters from the first "#" character onward are treated as
comment and discarded. See an exception in (6). Note that a
comment after a trailing "&" character will prevent the command from
continuing on the next line. Also note that for multi-line commands a
single leading "#" will comment out the entire command.
(3) The line is searched repeatedly for $ characters, which indicate
variables that are replaced with a text string. See an exception in
(6).
If the $ is followed by curly brackets, then the variable name is the
text inside the curly brackets. If no curly brackets follow the $,
then the variable name is the single character immediately following
the $. Thus $\{myTemp\} and $x refer to variable names "myTemp" and
"x".
If the $ is followed by parenthesis, then the text inside the
parenthesis is treated as an "immediate" variable and evaluated as an
"equal-style variable"_variable.html. This is a way to use numeric
formulas in an input script without having to assign them to variable
names. For example, these 3 input script lines:
variable X equal (xlo+xhi)/2+sqrt(v_area)
region 1 block $X 2 INF INF EDGE EDGE
variable X delete :pre
can be replaced by
region 1 block $((xlo+xhi)/2+sqrt(v_area)) 2 INF INF EDGE EDGE :pre
so that you do not have to define (or discard) a temporary variable X.
Note that neither the curly-bracket or immediate form of variables can
contain nested $ characters for other variables to substitute for.
Thus you cannot do this:
variable a equal 2
variable b2 equal 4
print "B2 = $\{b$a\}" :pre
Nor can you specify this $($x-1.0) for an immediate variable, but
you could use $(v_x-1.0), since the latter is valid syntax for an
"equal-style variable"_variable.html.
See the "variable"_variable.html command for more details of how
strings are assigned to variables and evaluated, and how they can be
used in input script commands.
(4) The line is broken into "words" separated by whitespace (tabs,
spaces). Note that words can thus contain letters, digits,
underscores, or punctuation characters.
(5) The first word is the command name. All successive words in the
line are arguments.
(6) If you want text with spaces to be treated as a single argument,
it can be enclosed in either double or single quotes. E.g.
print "Volume = $v"
print 'Volume = $v'
if "${steps} > 1000" then quit :pre
The quotes are removed when the single argument is stored internally.
See the "dump modify format"_dump_modify.html or "print"_print.html or
"if"_if.html commands for examples. A "#" or "$" character that is
between quotes will not be treated as a comment indicator in (2) or
substituted for as a variable in (3).
IMPORTANT NOTE: If the argument is itself a command that requires a
quoted argument (e.g. using a "print"_print.html command as part of an
"if"_if.html or "run every"_run.html command), then the double and
single quotes can be nested in the usual manner. See the doc pages
for those commands for examples. Only one of level of nesting is
allowed, but that should be sufficient for most use cases.
:line
3.3 Input script structure :h4,link(cmd_3)
This section describes the structure of a typical LAMMPS input script.
The "examples" directory in the LAMMPS distribution contains many
sample input scripts; the corresponding problems are discussed in
"Section_example"_Section_example.html, and animated on the "LAMMPS
WWW Site"_lws.
A LAMMPS input script typically has 4 parts:
Initialization
Atom definition
Settings
Run a simulation :ol
The last 2 parts can be repeated as many times as desired. I.e. run a
simulation, change some settings, run some more, etc. Each of the 4
parts is now described in more detail. Remember that almost all the
commands need only be used if a non-default value is desired.
(1) Initialization
Set parameters that need to be defined before atoms are created or
read-in from a file.
The relevant commands are "units"_units.html,
"dimension"_dimension.html, "newton"_newton.html,
"processors"_processors.html, "boundary"_boundary.html,
"atom_style"_atom_style.html, "atom_modify"_atom_modify.html.
If force-field parameters appear in the files that will be read, these
commands tell LAMMPS what kinds of force fields are being used:
"pair_style"_pair_style.html, "bond_style"_bond_style.html,
"angle_style"_angle_style.html, "dihedral_style"_dihedral_style.html,
"improper_style"_improper_style.html.
(2) Atom definition
There are 3 ways to define atoms in LAMMPS. Read them in from a data
or restart file via the "read_data"_read_data.html or
"read_restart"_read_restart.html commands. These files can contain
molecular topology information. Or create atoms on a lattice (with no
molecular topology), using these commands: "lattice"_lattice.html,
"region"_region.html, "create_box"_create_box.html,
"create_atoms"_create_atoms.html. The entire set of atoms can be
duplicated to make a larger simulation using the
"replicate"_replicate.html command.
(3) Settings
Once atoms and molecular topology are defined, a variety of settings
can be specified: force field coefficients, simulation parameters,
output options, etc.
Force field coefficients are set by these commands (they can also be
set in the read-in files): "pair_coeff"_pair_coeff.html,
"bond_coeff"_bond_coeff.html, "angle_coeff"_angle_coeff.html,
"dihedral_coeff"_dihedral_coeff.html,
"improper_coeff"_improper_coeff.html,
"kspace_style"_kspace_style.html, "dielectric"_dielectric.html,
"special_bonds"_special_bonds.html.
Various simulation parameters are set by these commands:
"neighbor"_neighbor.html, "neigh_modify"_neigh_modify.html,
"group"_group.html, "timestep"_timestep.html,
"reset_timestep"_reset_timestep.html, "run_style"_run_style.html,
"min_style"_min_style.html, "min_modify"_min_modify.html.
Fixes impose a variety of boundary conditions, time integration, and
diagnostic options. The "fix"_fix.html command comes in many flavors.
Various computations can be specified for execution during a
simulation using the "compute"_compute.html,
"compute_modify"_compute_modify.html, and "variable"_variable.html
commands.
Output options are set by the "thermo"_thermo.html, "dump"_dump.html,
and "restart"_restart.html commands.
(4) Run a simulation
A molecular dynamics simulation is run using the "run"_run.html
command. Energy minimization (molecular statics) is performed using
the "minimize"_minimize.html command. A parallel tempering
(replica-exchange) simulation can be run using the
"temper"_temper.html command.
:line
3.4 Commands listed by category :link(cmd_4),h4
This section lists all LAMMPS commands, grouped by category. The
"next section"_#cmd_5 lists the same commands alphabetically. Note
that some style options for some commands are part of specific LAMMPS
packages, which means they cannot be used unless the package was
included when LAMMPS was built. Not all packages are included in a
default LAMMPS build. These dependencies are listed as Restrictions
in the command's documentation.
Initialization:
"atom_modify"_atom_modify.html, "atom_style"_atom_style.html,
"boundary"_boundary.html, "dimension"_dimension.html,
"newton"_newton.html, "processors"_processors.html, "units"_units.html
Atom definition:
"create_atoms"_create_atoms.html, "create_box"_create_box.html,
"lattice"_lattice.html, "read_data"_read_data.html,
"read_dump"_read_dump.html, "read_restart"_read_restart.html,
"region"_region.html, "replicate"_replicate.html
Force fields:
"angle_coeff"_angle_coeff.html, "angle_style"_angle_style.html,
"bond_coeff"_bond_coeff.html, "bond_style"_bond_style.html,
"dielectric"_dielectric.html, "dihedral_coeff"_dihedral_coeff.html,
"dihedral_style"_dihedral_style.html,
"improper_coeff"_improper_coeff.html,
"improper_style"_improper_style.html,
"kspace_modify"_kspace_modify.html, "kspace_style"_kspace_style.html,
"pair_coeff"_pair_coeff.html, "pair_modify"_pair_modify.html,
"pair_style"_pair_style.html, "pair_write"_pair_write.html,
"special_bonds"_special_bonds.html
Settings:
"communicate"_communicate.html, "group"_group.html, "mass"_mass.html,
"min_modify"_min_modify.html, "min_style"_min_style.html,
"neigh_modify"_neigh_modify.html, "neighbor"_neighbor.html,
"reset_timestep"_reset_timestep.html, "run_style"_run_style.html,
"set"_set.html, "timestep"_timestep.html, "velocity"_velocity.html
Fixes:
"fix"_fix.html, "fix_modify"_fix_modify.html, "unfix"_unfix.html
Computes:
"compute"_compute.html, "compute_modify"_compute_modify.html,
"uncompute"_uncompute.html
Output:
"dump"_dump.html, "dump image"_dump_image.html,
"dump_modify"_dump_modify.html, "restart"_restart.html,
"thermo"_thermo.html, "thermo_modify"_thermo_modify.html,
"thermo_style"_thermo_style.html, "undump"_undump.html,
"write_data"_write_data.html, "write_restart"_write_restart.html
Actions:
"delete_atoms"_delete_atoms.html, "delete_bonds"_delete_bonds.html,
"displace_atoms"_displace_atoms.html, "change_box"_change_box.html,
"minimize"_minimize.html, "neb"_neb.html "prd"_prd.html,
"rerun"_rerun.html, "run"_run.html, "temper"_temper.html
Miscellaneous:
"clear"_clear.html, "echo"_echo.html, "if"_if.html,
"include"_include.html, "jump"_jump.html, "label"_label.html,
"log"_log.html, "next"_next.html, "print"_print.html,
"shell"_shell.html, "variable"_variable.html
:line
3.5 Individual commands :h4,link(cmd_5),link(comm)
This section lists all LAMMPS commands alphabetically, with a separate
listing below of styles within certain commands. The "previous
section"_#cmd_4 lists the same commands, grouped by category. Note
that some style options for some commands are part of specific LAMMPS
packages, which means they cannot be used unless the package was
included when LAMMPS was built. Not all packages are included in a
default LAMMPS build. These dependencies are listed as Restrictions
in the command's documentation.
"angle_coeff"_angle_coeff.html,
"angle_style"_angle_style.html,
"atom_modify"_atom_modify.html,
"atom_style"_atom_style.html,
"balance"_balance.html,
"bond_coeff"_bond_coeff.html,
"bond_style"_bond_style.html,
"boundary"_boundary.html,
"box"_box.html,
"change_box"_change_box.html,
"clear"_clear.html,
"communicate"_communicate.html,
"compute"_compute.html,
"compute_modify"_compute_modify.html,
"create_atoms"_create_atoms.html,
"create_box"_create_box.html,
"delete_atoms"_delete_atoms.html,
"delete_bonds"_delete_bonds.html,
"dielectric"_dielectric.html,
"dihedral_coeff"_dihedral_coeff.html,
"dihedral_style"_dihedral_style.html,
"dimension"_dimension.html,
"displace_atoms"_displace_atoms.html,
"dump"_dump.html,
"dump image"_dump_image.html,
"dump_modify"_dump_modify.html,
"echo"_echo.html,
"fix"_fix.html,
"fix_modify"_fix_modify.html,
"group"_group.html,
"if"_if.html,
"improper_coeff"_improper_coeff.html,
"improper_style"_improper_style.html,
"include"_include.html,
"jump"_jump.html,
"kspace_modify"_kspace_modify.html,
"kspace_style"_kspace_style.html,
"label"_label.html,
"lattice"_lattice.html,
"log"_log.html,
"mass"_mass.html,
"minimize"_minimize.html,
"min_modify"_min_modify.html,
"min_style"_min_style.html,
"neb"_neb.html,
"neigh_modify"_neigh_modify.html,
"neighbor"_neighbor.html,
"newton"_newton.html,
"next"_next.html,
"package"_package.html,
"pair_coeff"_pair_coeff.html,
"pair_modify"_pair_modify.html,
"pair_style"_pair_style.html,
"pair_write"_pair_write.html,
"partition"_partition.html,
"prd"_prd.html,
"print"_print.html,
"processors"_processors.html,
"quit"_quit.html,
"read_data"_read_data.html,
"read_dump"_read_dump.html,
"read_restart"_read_restart.html,
"region"_region.html,
"replicate"_replicate.html,
"rerun"_rerun.html,
"reset_timestep"_reset_timestep.html,
"restart"_restart.html,
"run"_run.html,
"run_style"_run_style.html,
"set"_set.html,
"shell"_shell.html,
"special_bonds"_special_bonds.html,
"suffix"_suffix.html,
"tad"_tad.html,
"temper"_temper.html,
"thermo"_thermo.html,
"thermo_modify"_thermo_modify.html,
"thermo_style"_thermo_style.html,
"timestep"_timestep.html,
"uncompute"_uncompute.html,
"undump"_undump.html,
"unfix"_unfix.html,
"units"_units.html,
"variable"_variable.html,
"velocity"_velocity.html,
"write_data"_write_data.html,
"write_restart"_write_restart.html :tb(c=6,ea=c)
These are commands contributed by users, which can be used if "LAMMPS
is built with the appropriate package"_Section_start.html#start_3.
"group2ndx"_group2ndx.html :tb(c=1,ea=c)
:line
Fix styles :h4
See the "fix"_fix.html command for one-line descriptions
of each style or click on the style itself for a full description:
"adapt"_fix_adapt.html,
"addforce"_fix_addforce.html,
"append/atoms"_fix_append_atoms.html,
"aveforce"_fix_aveforce.html,
"ave/atom"_fix_ave_atom.html,
"ave/correlate"_fix_ave_correlate.html,
"ave/histo"_fix_ave_histo.html,
"ave/spatial"_fix_ave_spatial.html,
"ave/time"_fix_ave_time.html,
"balance"_fix_balance.html,
"bond/break"_fix_bond_break.html,
"bond/create"_fix_bond_create.html,
"bond/swap"_fix_bond_swap.html,
"box/relax"_fix_box_relax.html,
"deform"_fix_deform.html,
"deposit"_fix_deposit.html,
"drag"_fix_drag.html,
"dt/reset"_fix_dt_reset.html,
"efield"_fix_efield.html,
"enforce2d"_fix_enforce2d.html,
"evaporate"_fix_evaporate.html,
"external"_fix_external.html,
"freeze"_fix_freeze.html,
"gcmc"_fix_gcmc.html,
"gld"_fix_gld.html,
"gravity"_fix_gravity.html,
"heat"_fix_heat.html,
"indent"_fix_indent.html,
"langevin"_fix_langevin.html,
"lineforce"_fix_lineforce.html,
"momentum"_fix_momentum.html,
"move"_fix_move.html,
"msst"_fix_msst.html,
"neb"_fix_neb.html,
"nph"_fix_nh.html,
"nphug"_fix_nphug.html,
"nph/asphere"_fix_nph_asphere.html,
"nph/sphere"_fix_nph_sphere.html,
"npt"_fix_nh.html,
"npt/asphere"_fix_npt_asphere.html,
"npt/sphere"_fix_npt_sphere.html,
"nve"_fix_nve.html,
"nve/asphere"_fix_nve_asphere.html,
"nve/asphere/noforce"_fix_nve_asphere_noforce.html,
"nve/body"_fix_nve_body.html,
"nve/limit"_fix_nve_limit.html,
"nve/line"_fix_nve_line.html,
"nve/noforce"_fix_nve_noforce.html,
"nve/sphere"_fix_nve_sphere.html,
"nve/tri"_fix_nve_tri.html,
"nvt"_fix_nh.html,
"nvt/asphere"_fix_nvt_asphere.html,
"nvt/sllod"_fix_nvt_sllod.html,
"nvt/sphere"_fix_nvt_sphere.html,
"orient/fcc"_fix_orient_fcc.html,
"planeforce"_fix_planeforce.html,
"poems"_fix_poems.html,
"pour"_fix_pour.html,
"press/berendsen"_fix_press_berendsen.html,
"print"_fix_print.html,
"property/atom"_fix_property_atom.html,
"qeq/comb"_fix_qeq_comb.html,
"reax/bonds"_fix_reax_bonds.html,
"recenter"_fix_recenter.html,
"restrain"_fix_restrain.html,
"rigid"_fix_rigid.html,
"rigid/nph"_fix_rigid.html,
"rigid/npt"_fix_rigid.html,
"rigid/nve"_fix_rigid.html,
"rigid/nvt"_fix_rigid.html,
"rigid/small"_fix_rigid.html,
"setforce"_fix_setforce.html,
"shake"_fix_shake.html,
"spring"_fix_spring.html,
"spring/rg"_fix_spring_rg.html,
"spring/self"_fix_spring_self.html,
"srd"_fix_srd.html,
"store/force"_fix_store_force.html,
"store/state"_fix_store_state.html,
"temp/berendsen"_fix_temp_berendsen.html,
"temp/rescale"_fix_temp_rescale.html,
"thermal/conductivity"_fix_thermal_conductivity.html,
"tmd"_fix_tmd.html,
"ttm"_fix_ttm.html,
"tune/kspace"_fix_tune_kspace.html,
"viscosity"_fix_viscosity.html,
"viscous"_fix_viscous.html,
"wall/colloid"_fix_wall.html,
"wall/gran"_fix_wall_gran.html,
"wall/harmonic"_fix_wall.html,
"wall/lj1043"_fix_wall.html,
"wall/lj126"_fix_wall.html,
"wall/lj93"_fix_wall.html,
"wall/piston"_fix_wall_piston.html,
"wall/reflect"_fix_wall_reflect.html,
"wall/region"_fix_wall_region.html,
"wall/srd"_fix_wall_srd.html :tb(c=8,ea=c)
These are fix styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"addtorque"_fix_addtorque.html,
"atc"_fix_atc.html,
"colvars"_fix_colvars.html,
"imd"_fix_imd.html,
"langevin/eff"_fix_langevin_eff.html,
"meso"_fix_meso.html,
"meso/stationary"_fix_meso_stationary.html,
"nph/eff"_fix_nh_eff.html,
"npt/eff"_fix_nh_eff.html,
"nve/eff"_fix_nve_eff.html,
"nvt/eff"_fix_nh_eff.html,
"nvt/sllod/eff"_fix_nvt_sllod_eff.html,
"phonon"_fix_phonon.html,
"qeq/reax"_fix_qeq_reax.html,
"reax/c/bonds"_fix_reax_bonds.html,
"reax/c/species"_fix_reaxc_species.html,
"smd"_fix_smd.html,
"temp/rescale/eff"_fix_temp_rescale_eff.html :tb(c=6,ea=c)
These are accelerated fix styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"freeze/cuda"_fix_freeze.html,
"addforce/cuda"_fix_addforce.html,
"aveforce/cuda"_fix_aveforce.html,
"enforce2d/cuda"_fix_enforce2d.html,
"gravity/cuda"_fix_gravity.html,
"gravity/omp"_fix_gravity.html,
"nph/omp"_fix_nh.html,
"nphug/omp"_fix_nphug.html,
"nph/asphere/omp"_fix_nph_asphere.html,
"nph/sphere/omp"_fix_nph_sphere.html,
"npt/cuda"_fix_nh.html,
"npt/omp"_fix_nh.html,
"npt/asphere/omp"_fix_npt_asphere.html,
"npt/sphere/omp"_fix_npt_sphere.html,
"nve/cuda"_fix_nh.html,
"nve/omp"_fix_nve.html,
"nve/sphere/omp"_fix_nve_sphere.html,
"nvt/cuda"_fix_nh.html,
"nvt/omp"_fix_nh.html,
"nvt/asphere/omp"_fix_nvt_asphere.html,
"nvt/sllod/omp"_fix_nvt_sllod.html,
"nvt/sphere/omp"_fix_nvt_sphere.html,
"qeq/comb/omp"_fix_qeq_comb.html,
"rigid/omp"_fix_rigid.html,
"rigid/nph/omp"_fix_rigid.html,
"rigid/npt/omp"_fix_rigid.html,
"rigid/nve/omp"_fix_rigid.html,
"rigid/nvt/omp"_fix_rigid.html,
"rigid/small/omp"_fix_rigid.html,
"setforce/cuda"_fix_setforce.html,
"shake/cuda"_fix_shake.html,
"temp/berendsen/cuda"_fix_temp_berendsen.html,
"temp/rescale/cuda"_fix_temp_rescale.html,
"temp/rescale/limit/cuda"_fix_temp_rescale.html,
"viscous/cuda"_fix_viscous.html :tb(c=6,ea=c)
:line
Compute styles :h4
See the "compute"_compute.html command for one-line descriptions of
each style or click on the style itself for a full description:
"angle/local"_compute_angle_local.html,
"atom/molecule"_compute_atom_molecule.html,
"body/local"_compute_body_local.html,
"bond/local"_compute_bond_local.html,
"centro/atom"_compute_centro_atom.html,
"cluster/atom"_compute_cluster_atom.html,
"cna/atom"_compute_cna_atom.html,
"com"_compute_com.html,
"com/molecule"_compute_com_molecule.html,
"contact/atom"_compute_contact_atom.html,
"coord/atom"_compute_coord_atom.html,
"damage/atom"_compute_damage_atom.html,
"dihedral/local"_compute_dihedral_local.html,
"displace/atom"_compute_displace_atom.html,
"erotate/asphere"_compute_erotate_asphere.html,
"erotate/rigid"_compute_erotate_rigid.html,
"erotate/sphere"_compute_erotate_sphere.html,
"erotate/sphere/atom"_compute_erotate_sphere_atom.html,
"event/displace"_compute_event_displace.html,
"group/group"_compute_group_group.html,
"gyration"_compute_gyration.html,
"gyration/molecule"_compute_gyration_molecule.html,
"heat/flux"_compute_heat_flux.html,
"improper/local"_compute_improper_local.html,
"inertia/molecule"_compute_inertia_molecule.html,
"ke"_compute_ke.html,
"ke/atom"_compute_ke_atom.html,
"ke/rigid"_compute_ke_rigid.html,
"msd"_compute_msd.html,
"msd/molecule"_compute_msd_molecule.html,
"msd/nongauss"_compute_msd_nongauss.html,
"pair"_compute_pair.html,
"pair/local"_compute_pair_local.html,
"pe"_compute_pe.html,
"pe/atom"_compute_pe_atom.html,
"pressure"_compute_pressure.html,
"property/atom"_compute_property_atom.html,
"property/local"_compute_property_local.html,
"property/molecule"_compute_property_molecule.html,
"rdf"_compute_rdf.html,
"reduce"_compute_reduce.html,
"reduce/region"_compute_reduce.html,
"slice"_compute_slice.html,
"stress/atom"_compute_stress_atom.html,
"temp"_compute_temp.html,
"temp/asphere"_compute_temp_asphere.html,
"temp/com"_compute_temp_com.html,
"temp/deform"_compute_temp_deform.html,
"temp/partial"_compute_temp_partial.html,
"temp/profile"_compute_temp_profile.html,
"temp/ramp"_compute_temp_ramp.html,
"temp/region"_compute_temp_region.html,
"temp/sphere"_compute_temp_sphere.html,
"ti"_compute_ti.html,
"voronoi/atom"_compute_voronoi_atom.html :tb(c=6,ea=c)
These are compute styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"ackland/atom"_compute_ackland_atom.html,
"basal/atom"_compute_basal_atom.html,
"ke/eff"_compute_ke_eff.html,
"ke/atom/eff"_compute_ke_atom_eff.html,
"meso_e/atom"_compute_meso_e_atom.html,
"meso_rho/atom"_compute_meso_rho_atom.html,
"meso_t/atom"_compute_meso_t_atom.html,
"temp/eff"_compute_temp_eff.html,
"temp/deform/eff"_compute_temp_deform_eff.html,
"temp/region/eff"_compute_temp_region_eff.html,
"temp/rotate"_compute_temp_rotate.html :tb(c=6,ea=c)
These are accelerated compute styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"pe/cuda"_compute_pe.html,
"pressure/cuda"_compute_pressure.html,
"temp/cuda"_compute_temp.html,
"temp/partial/cuda"_compute_temp_partial.html :tb(c=6,ea=c)
:line
Pair_style potentials :h4
See the "pair_style"_pair_style.html command for an overview of pair
potentials. Click on the style itself for a full description:
"none"_pair_none.html,
"hybrid"_pair_hybrid.html,
"hybrid/overlay"_pair_hybrid.html,
"adp"_pair_adp.html,
"airebo"_pair_airebo.html,
"beck"_pair_beck.html,
"body"_pair_body.html,
"bop"_pair_bop.html,
"born"_pair_born.html,
"born/coul/long"_pair_born.html,
"born/coul/msm"_pair_born.html,
"born/coul/wolf"_pair_born.html,
"brownian"_pair_brownian.html,
"brownian/poly"_pair_brownian.html,
"buck"_pair_buck.html,
"buck/coul/cut"_pair_buck.html,
"buck/coul/long"_pair_buck.html,
"buck/coul/msm"_pair_buck.html,
"buck/long/coul/long"_pair_buck_long.html,
"colloid"_pair_colloid.html,
"comb"_pair_comb.html,
"coul/cut"_pair_coul.html,
"coul/debye"_pair_coul.html,
"coul/dsf"_pair_coul.html,
"coul/long"_pair_coul.html,
"coul/msm"_pair_coul.html,
"coul/wolf"_pair_coul.html,
"dpd"_pair_dpd.html,
"dpd/tstat"_pair_dpd.html,
"dsmc"_pair_dsmc.html,
"eam"_pair_eam.html,
"eam/alloy"_pair_eam.html,
"eam/fs"_pair_eam.html,
"eim"_pair_eim.html,
"gauss"_pair_gauss.html,
"gayberne"_pair_gayberne.html,
"gran/hertz/history"_pair_gran.html,
"gran/hooke"_pair_gran.html,
"gran/hooke/history"_pair_gran.html,
"hbond/dreiding/lj"_pair_hbond_dreiding.html,
"hbond/dreiding/morse"_pair_hbond_dreiding.html,
"kim"_pair_kim.html,
"lcbop"_pair_lcbop.html,
"line/lj"_pair_line_lj.html,
"lj/charmm/coul/charmm"_pair_charmm.html,
"lj/charmm/coul/charmm/implicit"_pair_charmm.html,
"lj/charmm/coul/long"_pair_charmm.html,
"lj/charmm/coul/msm"_pair_charmm.html,
"lj/class2"_pair_class2.html,
"lj/class2/coul/cut"_pair_class2.html,
"lj/class2/coul/long"_pair_class2.html,
"lj/cut"_pair_lj.html,
"lj/cut/coul/cut"_pair_lj.html,
"lj/cut/coul/debye"_pair_lj.html,
"lj/cut/coul/dsf"_pair_lj.html,
"lj/cut/coul/long"_pair_lj.html,
"lj/cut/coul/msm"_pair_lj.html,
"lj/cut/dipole/cut"_pair_dipole.html,
"lj/cut/dipole/long"_pair_dipole.html,
"lj/cut/tip4p/cut"_pair_lj.html,
"lj/cut/tip4p/long"_pair_lj.html,
"lj/expand"_pair_lj_expand.html,
"lj/gromacs"_pair_gromacs.html,
"lj/gromacs/coul/gromacs"_pair_gromacs.html,
"lj/long/coul/long"_pair_lj_long.html,
"lj/long/dipole/long"_pair_dipole.html,
"lj/long/tip4p/long"_pair_lj_long.html,
"lj/smooth"_pair_lj_smooth.html,
"lj/smooth/linear"_pair_lj_smooth_linear.html,
"lj96/cut"_pair_lj96.html,
"lubricate"_pair_lubricate.html,
"lubricate/poly"_pair_lubricate.html,
"lubricateU"_pair_lubricateU.html,
"lubricateU/poly"_pair_lubricateU.html,
"meam"_pair_meam.html,
"mie/cut"_pair_mie.html,
"morse"_pair_morse.html,
"nm/cut"_pair_nm.html,
"nm/cut/coul/cut"_pair_nm.html,
"nm/cut/coul/long"_pair_nm.html,
"peri/lps"_pair_peri.html,
"peri/pmb"_pair_peri.html,
"peri/ves"_pair_peri.html,
"reax"_pair_reax.html,
"rebo"_pair_airebo.html,
"resquared"_pair_resquared.html,
"soft"_pair_soft.html,
"sw"_pair_sw.html,
"table"_pair_table.html,
"tersoff"_pair_tersoff.html,
"tersoff/mod"_pair_tersoff_mod.html,
"tersoff/zbl"_pair_tersoff_zbl.html,
"tip4p/cut"_pair_coul.html,
"tip4p/long"_pair_coul.html,
"tri/lj"_pair_tri_lj.html,
"yukawa"_pair_yukawa.html,
"yukawa/colloid"_pair_yukawa_colloid.html,
"zbl"_pair_zbl.html :tb(c=4,ea=c)
These are pair styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"awpmd/cut"_pair_awpmd.html,
"coul/diel"_pair_coul_diel.html,
"eam/cd"_pair_eam.html,
"edip"_pair_edip.html,
"eff/cut"_pair_eff.html,
"gauss/cut"_pair_gauss.html,
"list"_pair_list.html,
"lj/cut/dipole/sf"_pair_dipole.html,
"lj/sdk"_pair_sdk.html,
"lj/sdk/coul/long"_pair_sdk.html,
"lj/sdk/coul/msm"_pair_sdk.html,
"lj/sf"_pair_lj_sf.html,
"meam/spline"_pair_meam_spline.html,
"meam/sw/spline"_pair_meam_sw_spline.html,
"nb3b/harmonic"_pair_nb3b_harmonic.html,
"reax/c"_pair_reax_c.html,
"sph/heatconduction"_pair_sph_heatconduction.html,
"sph/idealgas"_pair_sph_idealgas.html,
"sph/lj"_pair_sph_lj.html,
"sph/rhosum"_pair_sph_rhosum.html,
"sph/taitwater"_pair_sph_taitwater.html,
"sph/taitwater/morris"_pair_sph_taitwater_morris.html,
"tersoff/table"_pair_tersoff.html :tb(c=4,ea=c)
These are accelerated pair styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"adp/omp"_pair_adp.html,
"airebo/omp"_pair_airebo.html,
"beck/gpu"_pair_beck.html,
"beck/omp"_pair_beck.html,
"born/coul/long/cuda"_pair_born.html,
"born/coul/long/gpu"_pair_born.html,
"born/coul/long/omp"_pair_born.html,
"born/coul/msm/omp"_pair_born.html,
"born/coul/wolf/gpu"_pair_born.html,
"born/coul/wolf/omp"_pair_born.html,
"born/gpu"_pair_born.html,
"born/omp"_pair_born.html,
"brownian/omp"_pair_brownian.html,
"brownian/poly/omp"_pair_brownian.html,
"buck/coul/cut/cuda"_pair_buck.html,
"buck/coul/cut/gpu"_pair_buck.html,
"buck/coul/cut/omp"_pair_buck.html,
"buck/coul/long/cuda"_pair_buck.html,
"buck/coul/long/gpu"_pair_buck.html,
"buck/coul/long/omp"_pair_buck.html,
"buck/coul/msm/omp"_pair_buck.html,
"buck/cuda"_pair_buck.html,
"buck/long/coul/long/omp"_pair_buck_long.html,
"buck/gpu"_pair_buck.html,
"buck/omp"_pair_buck.html,
"colloid/gpu"_pair_colloid.html,
"colloid/omp"_pair_colloid.html,
"comb/omp"_pair_comb.html,
"coul/cut/omp"_pair_coul.html,
"coul/debye/omp"_pair_coul.html,
"coul/dsf/gpu"_pair_coul.html,
"coul/long/gpu"_pair_coul.html,
"coul/long/omp"_pair_coul.html,
"coul/msm/omp"_pair_coul.html,
"coul/wolf"_pair_coul.html,
"dpd/omp"_pair_dpd.html,
"dpd/tstat/omp"_pair_dpd.html,
"eam/alloy/cuda"_pair_eam.html,
"eam/alloy/gpu"_pair_eam.html,
"eam/alloy/omp"_pair_eam.html,
"eam/alloy/opt"_pair_eam.html,
"eam/cd/omp"_pair_eam.html,
"eam/cuda"_pair_eam.html,
"eam/fs/cuda"_pair_eam.html,
"eam/fs/gpu"_pair_eam.html,
"eam/fs/omp"_pair_eam.html,
"eam/fs/opt"_pair_eam.html,
"eam/gpu"_pair_eam.html,
"eam/omp"_pair_eam.html,
"eam/opt"_pair_eam.html,
"edip/omp"_pair_edip.html,
"eim/omp"_pair_eim.html,
"gauss/gpu"_pair_gauss.html,
"gauss/omp"_pair_gauss.html,
"gayberne/gpu"_pair_gayberne.html,
"gayberne/omp"_pair_gayberne.html,
"gran/hertz/history/omp"_pair_gran.html,
"gran/hooke/cuda"_pair_gran.html,
"gran/hooke/history/omp"_pair_gran.html,
"gran/hooke/omp"_pair_gran.html,
"hbond/dreiding/lj/omp"_pair_hbond_dreiding.html,
"hbond/dreiding/morse/omp"_pair_hbond_dreiding.html,
"line/lj/omp"_pair_line_lj.html,
"lj/charmm/coul/charmm/cuda"_pair_charmm.html,
"lj/charmm/coul/charmm/omp"_pair_charmm.html,
"lj/charmm/coul/charmm/implicit/cuda"_pair_charmm.html,
"lj/charmm/coul/charmm/implicit/omp"_pair_charmm.html,
"lj/charmm/coul/long/cuda"_pair_charmm.html,
"lj/charmm/coul/long/gpu"_pair_charmm.html,
"lj/charmm/coul/long/omp"_pair_charmm.html,
"lj/charmm/coul/long/opt"_pair_charmm.html,
"lj/class2/coul/cut/cuda"_pair_class2.html,
"lj/class2/coul/cut/omp"_pair_class2.html,
"lj/class2/coul/long/cuda"_pair_class2.html,
"lj/class2/coul/long/gpu"_pair_class2.html,
"lj/class2/coul/long/omp"_pair_class2.html,
"lj/class2/coul/msm/omp"_pair_class2.html,
"lj/class2/cuda"_pair_class2.html,
"lj/class2/gpu"_pair_class2.html,
"lj/class2/omp"_pair_class2.html,
"lj/long/coul/long/omp"_pair_lj_long.html,
"lj/cut/coul/cut/cuda"_pair_lj.html,
"lj/cut/coul/cut/gpu"_pair_lj.html,
"lj/cut/coul/cut/omp"_pair_lj.html,
"lj/cut/coul/debye/cuda"_pair_lj.html,
"lj/cut/coul/debye/gpu"_pair_lj.html,
"lj/cut/coul/debye/omp"_pair_lj.html,
"lj/cut/coul/dsf/gpu"_pair_lj.html,
"lj/cut/coul/long/cuda"_pair_lj.html,
"lj/cut/coul/long/gpu"_pair_lj.html,
"lj/cut/coul/long/omp"_pair_lj.html,
"lj/cut/coul/long/opt"_pair_lj.html,
"lj/cut/coul/msm/gpu"_pair_lj.html,
"lj/cut/coul/msm/opt"_pair_lj.html,
"lj/cut/cuda"_pair_lj.html,
"lj/cut/dipole/cut/gpu"_pair_dipole.html,
"lj/cut/dipole/cut/omp"_pair_dipole.html,
"lj/cut/dipole/sf/gpu"_pair_dipole.html,
"lj/cut/dipole/sf/omp"_pair_dipole.html,
"lj/cut/experimental/cuda"_pair_lj.html,
"lj/cut/gpu"_pair_lj.html,
"lj/cut/omp"_pair_lj.html,
"lj/cut/opt"_pair_lj.html,
"lj/cut/tip4p/cut/omp"_pair_lj.html,
"lj/cut/tip4p/long/omp"_pair_lj.html,
"lj/cut/tip4p/long/opt"_pair_lj.html,
"lj/expand/cuda"_pair_lj_expand.html,
"lj/expand/gpu"_pair_lj_expand.html,
"lj/expand/omp"_pair_lj_expand.html,
"lj/gromacs/coul/gromacs/cuda"_pair_gromacs.html,
"lj/gromacs/coul/gromacs/omp"_pair_gromacs.html,
"lj/gromacs/cuda"_pair_gromacs.html,
"lj/gromacs/omp"_pair_gromacs.html,
"lj/long/coul/long/opt"_pair_lj_long.html,
"lj/sdk/gpu"_pair_sdk.html,
"lj/sdk/omp"_pair_sdk.html,
"lj/sdk/coul/long/gpu"_pair_sdk.html,
"lj/sdk/coul/long/omp"_pair_sdk.html,
"lj/sdk/coul/msm/omp"_pair_sdk.html,
"lj/sf/omp"_pair_lj_sf.html,
"lj/smooth/cuda"_pair_lj_smooth.html,
"lj/smooth/omp"_pair_lj_smooth.html,
"lj/smooth/linear/omp"_pair_lj_smooth_linear.html,
"lj96/cut/cuda"_pair_lj96.html,
"lj96/cut/gpu"_pair_lj96.html,
"lj96/cut/omp"_pair_lj96.html,
"lubricate/omp"_pair_lubricate.html,
"lubricate/poly/omp"_pair_lubricate.html,
"meam/spline/omp"_pair_meam_spline.html,
"mie/cut/gpu"_pair_mie.html,
"morse/cuda"_pair_morse.html,
"morse/gpu"_pair_morse.html,
"morse/omp"_pair_morse.html,
"morse/opt"_pair_morse.html,
"nb3b/harmonic/omp"_pair_nb3b_harmonic.html,
+"nm/cut/omp"_pair_nm.html,
+"nm/cut/coul/cut/omp"_pair_nm.html,
+"nm/cut/coul/long/omp"_pair_nm.html,
"peri/lps/omp"_pair_peri.html,
"peri/pmb/omp"_pair_peri.html,
"rebo/omp"_pair_airebo.html,
"resquared/gpu"_pair_resquared.html,
"resquared/omp"_pair_resquared.html,
"soft/gpu"_pair_soft.html,
"soft/omp"_pair_soft.html,
"sw/cuda"_pair_sw.html,
"sw/gpu"_pair_sw.html,
"sw/omp"_pair_sw.html,
"table/gpu"_pair_table.html,
"table/omp"_pair_table.html,
"tersoff/cuda"_pair_tersoff.html,
"tersoff/omp"_pair_tersoff.html,
+"tersoff/mod/omp"_pair_tersoff_mod.html,
"tersoff/table/omp"_pair_tersoff.html,
"tersoff/zbl/omp"_pair_tersoff_zbl.html,
"tip4p/cut/omp"_pair_coul.html,
"tip4p/long/omp"_pair_coul.html,
"tri/lj/omp"_pair_tri_lj.html,
"yukawa/gpu"_pair_yukawa.html,
"yukawa/omp"_pair_yukawa.html,
"yukawa/colloid/gpu"_pair_yukawa_colloid.html,
-"yukawa/colloid/omp"_pair_yukawa_colloid.html :tb(c=4,ea=c)
+"yukawa/colloid/omp"_pair_yukawa_colloid.html,
+"zbl/omp"_pair_zbl.html :tb(c=4,ea=c)
:line
Bond_style potentials :h4
See the "bond_style"_bond_style.html command for an overview of bond
potentials. Click on the style itself for a full description:
"none"_bond_none.html,
"hybrid"_bond_hybrid.html,
"class2"_bond_class2.html,
"fene"_bond_fene.html,
"fene/expand"_bond_fene_expand.html,
"harmonic"_bond_harmonic.html,
"morse"_bond_morse.html,
"nonlinear"_bond_nonlinear.html,
"quartic"_bond_quartic.html,
"table"_bond_table.html :tb(c=4,ea=c,w=100)
These are bond styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"harmonic/shift"_bond_harmonic_shift.html,
"harmonic/shift/cut"_bond_harmonic_shift_cut.html :tb(c=4,ea=c)
These are accelerated bond styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"class2/omp"_bond_class2.html,
"fene/omp"_bond_fene.html,
"fene/expand/omp"_bond_fene_expand.html,
"harmonic/omp"_bond_harmonic.html,
"harmonic/shift/omp"_bond_harmonic_shift.html,
"harmonic/shift/cut/omp"_bond_harmonic_shift_cut.html,
"morse/omp"_bond_morse.html,
"nonlinear/omp"_bond_nonlinear.html,
"quartic/omp"_bond_quartic.html,
"table/omp"_bond_table.html :tb(c=4,ea=c,w=100)
:line
Angle_style potentials :h4
See the "angle_style"_angle_style.html command for an overview of
angle potentials. Click on the style itself for a full description:
"none"_angle_none.html,
"hybrid"_angle_hybrid.html,
"charmm"_angle_charmm.html,
"class2"_angle_class2.html,
"cosine"_angle_cosine.html,
"cosine/delta"_angle_cosine_delta.html,
"cosine/periodic"_angle_cosine_periodic.html,
"cosine/squared"_angle_cosine_squared.html,
"harmonic"_angle_harmonic.html,
"table"_angle_table.html :tb(c=4,ea=c,w=100)
These are angle styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"sdk"_angle_sdk.html,
"cosine/shift"_angle_cosine_shift.html,
"cosine/shift/exp"_angle_cosine_shift_exp.html,
"dipole"_angle_dipole.html,
"fourier"_angle_fourier.html,
"fourier/simple"_angle_fourier_simple.html,
"quartic"_angle_quartic.html :tb(c=4,ea=c)
These are accelerated angle styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"charmm/omp"_angle_charmm.html,
"class2/omp"_angle_class2.html,
"cosine/omp"_angle_cosine.html,
"cosine/delta/omp"_angle_cosine_delta.html,
"cosine/periodic/omp"_angle_cosine_periodic.html,
"cosine/shift/omp"_angle_cosine_shift.html,
"cosine/shift/exp/omp"_angle_cosine_shift_exp.html,
"cosine/squared/omp"_angle_cosine_squared.html,
"dipole/omp"_angle_dipole.html
"fourier/omp"_angle_fourier.html,
"fourier/simple/omp"_angle_fourier_simple.html,
"harmonic/omp"_angle_harmonic.html,
"quartic/omp"_angle_quartic.html
"table/omp"_angle_table.html :tb(c=4,ea=c,w=100)
:line
Dihedral_style potentials :h4
See the "dihedral_style"_dihedral_style.html command for an overview
of dihedral potentials. Click on the style itself for a full
description:
"none"_dihedral_none.html,
"hybrid"_dihedral_hybrid.html,
"charmm"_dihedral_charmm.html,
"class2"_dihedral_class2.html,
"harmonic"_dihedral_harmonic.html,
"helix"_dihedral_helix.html,
"multi/harmonic"_dihedral_multi_harmonic.html,
"opls"_dihedral_opls.html :tb(c=4,ea=c,w=100)
These are dihedral styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"cosine/shift/exp"_dihedral_cosine_shift_exp.html,
"fourier"_dihedral_fourier.html,
"nharmonic"_dihedral_nharmonic.html,
"quadratic"_dihedral_quadratic.html,
"table"_dihedral_table.html :tb(c=4,ea=c)
These are accelerated dihedral styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"charmm/omp"_dihedral_charmm.html,
"class2/omp"_dihedral_class2.html,
"cosine/shift/exp/omp"_dihedral_cosine_shift_exp.html,
"fourier/omp"_dihedral_fourier.html,
"harmonic/omp"_dihedral_harmonic.html,
"helix/omp"_dihedral_helix.html,
"multi/harmonic/omp"_dihedral_multi_harmonic.html,
"nharmonic/omp"_dihedral_nharmonic.html,
"opls/omp"_dihedral_opls.html
"quadratic/omp"_dihedral_quadratic.html,
"table/omp"_dihedral_table.html :tb(c=4,ea=c,w=100)
:line
Improper_style potentials :h4
See the "improper_style"_improper_style.html command for an overview
of improper potentials. Click on the style itself for a full
description:
"none"_improper_none.html,
"hybrid"_improper_hybrid.html,
"class2"_improper_class2.html,
"cvff"_improper_cvff.html,
"harmonic"_improper_harmonic.html,
"umbrella"_improper_umbrella.html :tb(c=4,ea=c,w=100)
These are improper styles contributed by users, which can be used if
"LAMMPS is built with the appropriate
package"_Section_start.html#start_3.
"cossq"_improper_cossq.html,
"fourier"_improper_fourier.html,
"ring"_improper_ring.html :tb(c=4,ea=c)
These are accelerated improper styles, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"class2/omp"_improper_class2.html,
"cossq/omp"_improper_cossq.html,
"cvff/omp"_improper_cvff.html,
"fourier/omp"_improper_fourier.html,
"harmonic/omp"_improper_harmonic.html,
"ring/omp"_improper_ring.html,
"umbrella/omp"_improper_umbrella.html :tb(c=4,ea=c,w=100)
:line
Kspace solvers :h4
See the "kspace_style"_kspace_style.html command for an overview of
Kspace solvers. Click on the style itself for a full description:
"ewald"_kspace_style.html,
"ewald/disp"_kspace_style.html,
"msm"_kspace_style.html,
"msm/cg"_kspace_style.html,
"pppm"_kspace_style.html,
"pppm/cg"_kspace_style.html,
"pppm/disp"_kspace_style.html,
"pppm/disp/tip4p"_kspace_style.html,
"pppm/tip4p"_kspace_style.html :tb(c=4,ea=c,w=100)
These are accelerated Kspace solvers, which can be used if LAMMPS is
built with the "appropriate accelerated
package"_Section_accelerate.html.
"ewald/omp"_kspace_style.html,
"msm/omp"_kspace_style.html,
"msm/cg/omp"_kspace_style.html,
"pppm/cuda"_kspace_style.html,
"pppm/gpu"_kspace_style.html,
"pppm/omp"_kspace_style.html,
"pppm/cg/omp"_kspace_style.html,
"pppm/tip4p/omp"_kspace_style.html :tb(c=4,ea=c)
diff --git a/doc/pair_tersoff_mod.txt b/doc/pair_tersoff_mod.txt
index b4b678a53..1b3bcd51d 100644
--- a/doc/pair_tersoff_mod.txt
+++ b/doc/pair_tersoff_mod.txt
@@ -1,167 +1,190 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style tersoff/mod command :h3
+pair_style tersoff/mod/omp command :h3
[Syntax:]
pair_style tersoff/mod :pre
[Examples:]
pair_style tersoff/mod
pair_coeff * * Si.tersoff.mod Si Si :pre
[Description:]
The {tersoff/mod} style computes a bond-order type interatomic
potential "(Kumagai)"_#Kumagai based on a 3-body Tersoff potential
"(Tersoff_1)"_#Tersoff_1, "(Tersoff_2)"_#Tersoff_2 with modified
cutoff function and angular-dependent term, giving the energy E of a
system of atoms as
:c,image(Eqs/pair_tersoff_mod.jpg)
where f_R is a two-body term and f_A includes three-body interactions.
The summations in the formula are over all neighbors J and K of atom I
within a cutoff distance = R + D.
The modified cutoff function f_C proposed by "(Murty)"_#Murty and
having a continuous second-order differential is employed. The
angular-dependent term g(theta) was modified to increase the
flexibility of the potential.
The {tersoff/mod} potential is fitted to both the elastic constants
and melting point by employing the modified Tersoff potential function
form in which the angular-dependent term is improved. The model
performs extremely well in describing the crystalline, liquid, and
amorphous phases "(Schelling)"_#Schelling.
Only a single pair_coeff command is used with the {tersoff/mod} style
which specifies a Tersoff/MOD potential file with parameters for all
needed elements. These are mapped to LAMMPS atom types by specifying
N additional arguments after the filename in the pair_coeff command,
where N is the number of LAMMPS atom types:
filename
N element names = mapping of Tersoff/MOD elements to atom types :ul
As an example, imagine the Si.tersoff_mod file has Tersoff values for Si.
If your LAMMPS simulation has 3 Si atoms types, you would use the following
pair_coeff command:
pair_coeff * * Si.tersoff_mod Si Si Si :pre
The 1st 2 arguments must be * * so as to span all LAMMPS atom types.
The three Si arguments map LAMMPS atom types 1,2,3 to the Si element
in the Tersoff/MOD file. If a mapping value is specified as NULL, the
mapping is not performed. This can be used when a {tersoff/mod}
potential is used as part of the {hybrid} pair style. The NULL values
are placeholders for atom types that will be used with other
potentials.
Tersoff/MOD file in the {potentials} directory of the LAMMPS
distribution have a ".tersoff.mod" suffix. Lines that are not blank
or comments (starting with #) define parameters for a triplet of
elements. The parameters in a single entry correspond to coefficients
in the formula above:
element 1 (the center atom in a 3-body interaction)
element 2 (the atom bonded to the center atom)
element 3 (the atom influencing the 1-2 bond in a bond-order sense)
beta
alpha
h
eta
beta_ters = 1 (dummy parameter)
lambda2 (1/distance units)
B (energy units)
R (distance units)
D (distance units)
lambda1 (1/distance units)
A (energy units)
n
c1
c2
c3
c4
c5 :ul
The n, eta, lambda2, B, lambda1, and A parameters are only used for
two-body interactions. The beta, alpha, c1, c2, c3, c4, c5, h
parameters are only used for three-body interactions. The R and D
parameters are used for both two-body and three-body interactions. The
non-annotated parameters are unitless.
The Tersoff/MOD potential file must contain entries for all the elements
listed in the pair_coeff command. It can also contain entries for
additional elements not being used in a particular simulation; LAMMPS
ignores those entries.
For a single-element simulation, only a single entry is required
(e.g. SiSiSi). As annotated above, the first element in the entry is
the center atom in a three-body interaction and it is bonded to the
2nd atom and the bond is influenced by the 3rd atom. Thus an entry
for SiSiSi means Si bonded to a Si with another Si atom influencing the bond.
:line
+Styles with a {cuda}, {gpu}, {omp}, or {opt} suffix are functionally
+the same as the corresponding style without the suffix. They have
+been optimized to run faster, depending on your available hardware, as
+discussed in "Section_accelerate"_Section_accelerate.html of the
+manual. The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the USER-CUDA, GPU, USER-OMP and OPT
+packages, respectively. They are only enabled if LAMMPS was built with
+those packages. See the "Making LAMMPS"_Section_start.html#start_3
+section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
+use the "suffix"_suffix.html command in your input script.
+
+See "Section_accelerate"_Section_accelerate.html of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+:line
+
[Mixing, shift, table, tail correction, restart, rRESPA info]:
This pair style does not support the "pair_modify"_pair_modify.html
shift, table, and tail options.
This pair style does not write its information to "binary restart
files"_restart.html, since it is stored in potential files. Thus, you
need to re-specify the pair_style and pair_coeff commands in an input
script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
:line
[Restrictions:]
This pair style is part of the MANYBODY package. It is only enabled
if LAMMPS was built with that package (which it is by default). See
the "Making LAMMPS"_Section_start.html#start_3 section for more info.
This pair style requires the "newton"_newton.html setting to be "on"
for pair interactions.
The Tersoff/MOD potential files provided with LAMMPS (see the potentials
directory) are parameterized for metal "units"_units.html. You can
use the Tersoff/MOD potential with any LAMMPS units, but you would need to
create your own Tersoff/MOD potential file with coefficients listed in the
appropriate units if your simulation doesn't use "metal" units.
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:line
:link(Kumagai)
[(Kumagai)] T. Kumagai, S. Izumi, S. Hara, S. Sakai,
Comp. Mat. Science, 39, 457 (2007).
:link(Tersoff_1)
[(Tersoff_1)] J. Tersoff, Phys Rev B, 37, 6991 (1988).
:link(Tersoff_2)
[(Tersoff_2)] J. Tersoff, Phys Rev B, 38, 9902 (1988).
:link(Murty)
[(Murty)] M.V.R. Murty, H.A. Atwater, Phys Rev B, 51, 4889 (1995).
:link(Schelling)
[(Schelling)] Patrick K. Schelling, Comp. Mat. Science, 44, 274 (2008).
diff --git a/doc/pair_zbl.txt b/doc/pair_zbl.txt
index 34a949a83..b9c74a492 100644
--- a/doc/pair_zbl.txt
+++ b/doc/pair_zbl.txt
@@ -1,106 +1,129 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style zbl command :h3
+pair_style zbl/omp command :h3
[Syntax:]
pair_style zbl inner outer :pre
inner = distance where switching function begins
outer = global cutoff for ZBL interaction :ul
[Examples:]
pair_style zbl 3.0 4.0
pair_coeff * * 73.0
pair_coeff 1 1 14.0 :pre
[Description:]
Style {zbl} computes the Ziegler-Biersack-Littmark (ZBL) screened nuclear
repulsion for describing high-energy collisions between atoms.
"(Ziegler)"_#Ziegler. It includes an additional switching function
that ramps the energy, force, and curvature smoothly to zero
between an inner and outer cutoff. The potential
energy due to a pair of atoms at a distance r_ij is given by:
:c,image(Eqs/pair_zbl.jpg)
where e is the electron
charge, epsilon_0 is the electrical permittivity of vacuum, and
Z_i and Z_j are the nuclear charges of the two atoms in electron
charge units.
The switching
function S(r) is identical to that used by
"pair_style lj/gromacs"_pair_gromacs.html.
Here, the inner and outer cutoff are the same
for all pairs of atom types.
The following coefficient must be defined for each pair of atom types
via the "pair_coeff"_pair_coeff.html command as in the examples above,
or in the LAMMPS data file.
Z can not be specified for two different atoms types.
Therefore the lists of atom types I and atom types J must match.
Z (electron charge) :ul
Although Z must be defined for all atom type pairs I,J, it is only
stored for individual atom types, i.e. when I = J. Z is normally equal
to the atomic number of the atom type.
IMPORTANT NOTE: The numerical values of the exponential decay
constants in the screening function depend on the unit of distance. In
the above equation they are given for units of angstroms. LAMMPS will
automatically convert these values to the distance unit of the
specified LAMMPS "units"_units.html setting. The values of Z should
always be given in units of electron charge.
:line
+Styles with a {cuda}, {gpu}, {omp}, or {opt} suffix are functionally
+the same as the corresponding style without the suffix. They have
+been optimized to run faster, depending on your available hardware, as
+discussed in "Section_accelerate"_Section_accelerate.html of the
+manual. The accelerated styles take the same arguments and should
+produce the same results, except for round-off and precision issues.
+
+These accelerated styles are part of the USER-CUDA, GPU, USER-OMP and OPT
+packages, respectively. They are only enabled if LAMMPS was built with
+those packages. See the "Making LAMMPS"_Section_start.html#start_3
+section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Section_start.html#start_7 when you invoke LAMMPS, or you can
+use the "suffix"_suffix.html command in your input script.
+
+See "Section_accelerate"_Section_accelerate.html of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+:line
+
[Mixing, shift, table, tail correction, restart, rRESPA info]:
Mixing is not relevant for this pair style, since as explained above,
Z values are stored on a per-type basis, and both Zi and Zj are used
explicitly in the ZBL formula.
The ZBL pair style does not support the "pair_modify"_pair_modify.html
shift option, since the ZBL interaction is already smoothed to 0.0 at
the cutoff.
The "pair_modify"_pair_modify.html table option is not relevant for
this pair style.
This pair style does not support the "pair_modify"_pair_modify.html
tail option for adding long-range tail corrections to energy and
pressure, since there are no corrections for a potential that goes to
0.0 at the cutoff.
This pair style does not write information to "binary restart
files"_restart.html, so pair_style and pair_coeff commands must be
specified in an input script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
:line
[Restrictions:] none
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:line
:link(Ziegler)
[(Ziegler)] J.F. Ziegler, J. P. Biersack and U. Littmark, "The
Stopping and Range of Ions in Matter," Volume 1, Pergamon, 1985.
diff --git a/doc/shell.txt b/doc/shell.txt
index 1369df80a..793a9fc8d 100644
--- a/doc/shell.txt
+++ b/doc/shell.txt
@@ -1,98 +1,106 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
shell command :h3
[Syntax:]
shell cmd args :pre
-cmd = {cd} or {mkdir} or {mv} or {rm} or {rmdir} or arbitrary command :ulb,l
+cmd = {cd} or {mkdir} or {mv} or {rm} or {rmdir} or {setenv} or arbitrary command :ulb,l
{cd} arg = dir
dir = directory to change to
{mkdir} args = dir1 dir2 ...
dir1,dir2 = one or more directories to create
{mv} args = old new
old = old filename
new = new filename
{rm} args = file1 file2 ...
file1,file2 = one or more filenames to delete
{rmdir} args = dir1 dir2 ...
dir1,dir2 = one or more directories to delete
+ {setenv} args = name value
+ name = name of the environment variable
+ value = value of the environment variable
anything else is passed as a command to the shell for direct execution :pre
:ule
[Examples:]
shell cd sub1
shell cd ..
shell mkdir tmp1 tmp2 tmp3
shell rmdir tmp1
shell mv log.lammps hold/log.1
shell rm TMP/file1 TMP/file2
+shell setenv LAMMPS_POTENTIALS ../../potentials
shell my_setup file1 10 file2
shell my_post_process 100 dump.out :pre
[Description:]
Execute a shell command. A few simple file-based shell commands are
supported directly, in Unix-style syntax. Any command not listed
above is passed as-is to the C-library system() call, which invokes
the command in a shell.
This is means to invoke other commands from your input script. For
example, you can move files around in preparation for the next section
of the input script. Or you can run a program that pre-processes data
for input into LAMMPS. Or you can run a program that post-processes
LAMMPS output data.
With the exception of {cd}, all commands, including ones invoked via a
system() call, are executed by only a single processor, so that
files/directories are not being manipulated by multiple processors.
The {cd} cmd executes the Unix "cd" command to change the working
directory. All subsequent LAMMPS commands that read/write files will
use the new directory. All processors execute this command.
The {mkdir} cmd executes the Unix "mkdir" command to create one or
more directories.
The {mv} cmd executes the Unix "mv" command to rename a file and/or
move it to a new directory.
The {rm} cmd executes the Unix "rm" command to remove one or more
files.
The {rmdir} cmd executes the Unix "rmdir" command to remove one or
more directories. A directory must be empty to be successfully
removed.
+The {setenv} cmd defines or updates an environment variable directly.
+Since this command does not pass through the shell, no shell variable
+expansion or globbing is performed, but the value string used literally.
+
Any other cmd is passed as-is to the shell along with its arguments as
one string, invoked by the C-library system() call. For example,
these lines in your input script:
variable n equal 10
variable foo string file2
shell my_setup file1 $n $\{foo\} :pre
would be the same as invoking
% my_setup file1 10 file2 :pre
from a command-line prompt. The executable program "my_setup" is run
with 3 arguments: file1 10 file2.
[Restrictions:]
LAMMPS does not detect errors or print warnings when any of these
commands execute. E.g. if the specified directory does not exist,
executing the {cd} command will silently do nothing.
[Related commands:] none
[Default:] none
diff --git a/src/MAKE/Makefile.openmpi-omp b/src/MAKE/Makefile.openmpi-omp
index 7dc125dc5..2928dbedf 100644
--- a/src/MAKE/Makefile.openmpi-omp
+++ b/src/MAKE/Makefile.openmpi-omp
@@ -1,105 +1,105 @@
# openmpi-omp = Fedora 12, gcc-4.4.3 /w OpenMP, OpenMPI-1.4.1, FFTW-2.1.5
# to select openmpi (over mpich) do: module load openmpi-x86_64
SHELL = /bin/sh
# ---------------------------------------------------------------------
# compiler/linker settings
# specify flags and libraries needed for your compiler
OPENMP= -fopenmp
CC = mpic++ $(OPENMP)
CCFLAGS = -O3 -fomit-frame-pointer -fno-rtti -fno-exceptions \
-march=native -ffast-math -mpc64 -g -fstrict-aliasing \
-Wall -W -Wno-unused-parameter -Wno-maybe-uninitialized
DEPFLAGS = -MM
LINK = mpic++ $(OPENMP)
LINKFLAGS = -O -g -fno-rtti -fno-exceptions -mpc64
LIB =
ARCHIVE = ar
ARFLAGS = -rcsv
SIZE = size
# ---------------------------------------------------------------------
# LAMMPS-specific settings
# specify settings for LAMMPS features you will use
# if you change any -D setting, do full re-compile after "make clean"
# LAMMPS ifdef settings, OPTIONAL
# see possible settings in doc/Section_start.html#2_2 (step 4)
-LMP_INC = -DLAMMPS_GZIP -DLAMMPS_MEMALIGN=32 \
+LMP_INC = -DLAMMPS_GZIP -DLAMMPS_MEMALIGN=64 \
-DLAMMPS_JPEG
# MPI library, REQUIRED
# see discussion in doc/Section_start.html#2_2 (step 5)
# can point to dummy MPI library in src/STUBS as in Makefile.serial
# INC = path for mpi.h, MPI compiler settings
# PATH = path for MPI library
# LIB = name of MPI library
MPI_INC = -DOMPI_SKIP_MPICXX=1
MPI_PATH =
MPI_LIB =
# FFT library, OPTIONAL
# see discussion in doc/Section_start.html#2_2 (step 6)
# can be left blank to use provided KISS FFT library
# INC = -DFFT setting, e.g. -DFFT_FFTW, FFT compiler settings
# PATH = path for FFT library
# LIB = name of FFT library
FFT_INC = -DFFT_FFTW3
FFT_PATH =
FFT_LIB = -lfftw3
# JPEG library, OPTIONAL
# see discussion in doc/Section_start.html#2_2 (step 7)
# only needed if -DLAMMPS_JPEG listed with LMP_INC
# INC = path for jpeglib.h
# PATH = path for JPEG library
# LIB = name of JPEG library
JPG_INC =
JPG_PATH =
JPG_LIB = -ljpeg
# ---------------------------------------------------------------------
# build rules and dependencies
# no need to edit this section
include Makefile.package.settings
include Makefile.package
EXTRA_INC = $(LMP_INC) $(PKG_INC) $(MPI_INC) $(FFT_INC) $(JPG_INC) $(PKG_SYSINC)
EXTRA_PATH = $(PKG_PATH) $(MPI_PATH) $(FFT_PATH) $(JPG_PATH) $(PKG_SYSPATH)
EXTRA_LIB = $(PKG_LIB) $(MPI_LIB) $(FFT_LIB) $(JPG_LIB) $(PKG_SYSLIB)
# Path to src files
vpath %.cpp ..
vpath %.h ..
# Link target
$(EXE): $(OBJ)
$(LINK) $(LINKFLAGS) $(EXTRA_PATH) $(OBJ) $(EXTRA_LIB) $(LIB) -o $(EXE)
$(SIZE) $(EXE)
# Library target
lib: $(OBJ)
$(ARCHIVE) $(ARFLAGS) $(EXE) $(OBJ)
# Compilation rules
%.o:%.cpp
$(CC) $(CCFLAGS) $(EXTRA_INC) -c $<
%.d:%.cpp
$(CC) $(CCFLAGS) $(EXTRA_INC) $(DEPFLAGS) $< > $@
# Individual dependencies
DEPENDS = $(OBJ:.o=.d)
sinclude $(DEPENDS)
diff --git a/src/MANYBODY/pair_tersoff_mod.cpp b/src/MANYBODY/pair_tersoff_mod.cpp
index f29b8c0c5..b452be0fa 100644
--- a/src/MANYBODY/pair_tersoff_mod.cpp
+++ b/src/MANYBODY/pair_tersoff_mod.cpp
@@ -1,361 +1,361 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) 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.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL) - original Tersoff implementation
Vitaly Dozhdikov (JIHT of RAS) - MOD addition
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_tersoff_mod.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairTersoffMOD::PairTersoffMOD(LAMMPS *lmp) : PairTersoff(lmp) {}
/* ---------------------------------------------------------------------- */
void PairTersoffMOD::read_file(char *file)
{
int params_per_line = 20;
char **words = new char*[params_per_line+1];
memory->sfree(params);
params = NULL;
nparams = maxparam = 0;
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
- fp = fopen(file,"r");
+ fp = open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open Tersoff potential file %s",file);
error->one(FLERR,str);
}
}
// read each line out of file, skipping blank lines or leading '#'
// store line of params if all 3 element tags are in element list
int n,nwords,ielement,jelement,kelement;
char line[MAXLINE],*ptr;
int eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if (ptr = strchr(line,'#')) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if (ptr = strchr(line,'#')) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in Tersoff potential file");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next line
for (ielement = 0; ielement < nelements; ielement++)
if (strcmp(words[0],elements[ielement]) == 0) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (strcmp(words[1],elements[jelement]) == 0) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (strcmp(words[2],elements[kelement]) == 0) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].powerm = atof(words[3]);
params[nparams].lam3 = atof(words[4]);
params[nparams].h = atof(words[5]);
params[nparams].powern = atof(words[6]);
params[nparams].beta = atof(words[7]);
params[nparams].lam2 = atof(words[8]);
params[nparams].bigb = atof(words[9]);
params[nparams].bigr = atof(words[10]);
params[nparams].bigd = atof(words[11]);
params[nparams].lam1 = atof(words[12]);
params[nparams].biga = atof(words[13]);
params[nparams].powern_del = atof(words[14]);
params[nparams].c1 = atof(words[15]);
params[nparams].c2 = atof(words[16]);
params[nparams].c3 = atof(words[17]);
params[nparams].c4 = atof(words[18]);
params[nparams].c5 = atof(words[19]);
// currently only allow m exponent of 1
params[nparams].powermint = int(params[nparams].powerm);
if (
params[nparams].lam3 < 0.0 || params[nparams].powern < 0.0 ||
params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||
params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
params[nparams].bigd < 0.0 ||
params[nparams].bigd > params[nparams].bigr ||
params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].powerm - params[nparams].powermint != 0.0 ||
(params[nparams].powermint != 3 && params[nparams].powermint != 1))
error->all(FLERR,"Illegal Tersoff parameter");
nparams++;
}
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairTersoffMOD::setup()
{
int i,j,k,m,n;
// set elem2param for all element triplet combinations
// must be a single exact match to lines read from file
// do not allow for ACB in place of ABC
memory->destroy(elem2param);
memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param");
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement &&
k == params[m].kelement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
elem2param[i][j][k] = n;
}
// compute parameter values derived from inputs
for (m = 0; m < nparams; m++) {
params[m].cut = params[m].bigr + params[m].bigd;
params[m].cutsq = params[m].cut*params[m].cut;
params[m].ca1 = pow(2.0*params[m].powern_del*1.0e-16,-1.0/params[m].powern);
params[m].ca4 = 1.0/params[m].ca1;
}
// set cutmax to max of all params
cutmax = 0.0;
for (m = 0; m < nparams; m++)
if (params[m].cut > cutmax) cutmax = params[m].cut;
}
/* ---------------------------------------------------------------------- */
double PairTersoffMOD::zeta(Param *param, double rsqij, double rsqik,
double *delrij, double *delrik)
{
double rij,rik,costheta,arg,ex_delr;
rij = sqrt(rsqij);
rik = sqrt(rsqik);
costheta = (delrij[0]*delrik[0] + delrij[1]*delrik[1] +
delrij[2]*delrik[2]) / (rij*rik);
if (param->powermint == 3) arg = pow(param->lam3 * (rij-rik),3.0);
else arg = param->lam3 * (rij-rik);
if (arg > 69.0776) ex_delr = 1.e30;
else if (arg < -69.0776) ex_delr = 0.0;
else ex_delr = exp(arg);
return ters_fc(rik,param) * ters_gijk_mod(costheta,param) * ex_delr;
}
/* ---------------------------------------------------------------------- */
double PairTersoffMOD::ters_fc(double r, Param *param)
{
double ters_R = param->bigr;
double ters_D = param->bigd;
if (r < ters_R-ters_D) return 1.0;
if (r > ters_R+ters_D) return 0.0;
return 0.5*(1.0 - 1.125*sin(MY_PI2*(r - ters_R)/ters_D) -
0.125*sin(3*MY_PI2*(r - ters_R)/ters_D));
}
/* ---------------------------------------------------------------------- */
double PairTersoffMOD::ters_fc_d(double r, Param *param)
{
double ters_R = param->bigr;
double ters_D = param->bigd;
if (r < ters_R-ters_D) return 0.0;
if (r > ters_R+ters_D) return 0.0;
return -(0.375*MY_PI4/ters_D) * (3*cos(MY_PI2*(r - ters_R)/ters_D) +
cos(3*MY_PI2*(r - ters_R)/ters_D));
}
/* ---------------------------------------------------------------------- */
double PairTersoffMOD::ters_bij(double zeta, Param *param)
{
double tmp = param->beta * zeta;
if (tmp > param->ca1) return pow(tmp, -param->powern/(2.0*param->powern_del));
if (tmp < param->ca4) return 1.0;
return pow(1.0 + pow(tmp,param->powern), -1.0/(2.0*param->powern_del));
}
/* ---------------------------------------------------------------------- */
double PairTersoffMOD::ters_bij_d(double zeta, Param *param)
{
double tmp = param->beta * zeta;
if (tmp > param->ca1) return -0.5*(param->powern/param->powern_del)*
pow(tmp,-0.5*(param->powern/param->powern_del)) / zeta;
if (tmp < param->ca4) return 0.0;
double tmp_n = pow(tmp,param->powern);
return -0.5 *(param->powern/param->powern_del)*
pow(1.0+tmp_n, -1.0-(1.0/(2.0*param->powern_del)))*tmp_n / zeta;
}
/* ---------------------------------------------------------------------- */
void PairTersoffMOD::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk,
Param *param)
{
double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
double dcosdri[3],dcosdrj[3],dcosdrk[3];
fc = ters_fc(rik,param);
dfc = ters_fc_d(rik,param);
if (param->powermint == 3) tmp = pow(param->lam3 * (rij-rik),3.0);
else tmp = param->lam3 * (rij-rik);
if (tmp > 69.0776) ex_delr = 1.e30;
else if (tmp < -69.0776) ex_delr = 0.0;
else ex_delr = exp(tmp);
if (param->powermint == 3)
ex_delr_d = 3.0*pow(param->lam3,3.0) * pow(rij-rik,2.0)*ex_delr;
else ex_delr_d = param->lam3 * ex_delr;
cos_theta = vec3_dot(rij_hat,rik_hat);
gijk = ters_gijk_mod(cos_theta,param);
gijk_d = ters_gijk_d_mod(cos_theta,param);
costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk);
// compute the derivative wrt Ri
// dri = -dfc*gijk*ex_delr*rik_hat;
// dri += fc*gijk_d*ex_delr*dcosdri;
// dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat);
vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri);
vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri);
vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri);
vec3_scale(prefactor,dri,dri);
// compute the derivative wrt Rj
// drj = fc*gijk_d*ex_delr*dcosdrj;
// drj += fc*gijk*ex_delr_d*rij_hat;
vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj);
vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj);
vec3_scale(prefactor,drj,drj);
// compute the derivative wrt Rk
// drk = dfc*gijk*ex_delr*rik_hat;
// drk += fc*gijk_d*ex_delr*dcosdrk;
// drk += -fc*gijk*ex_delr_d*rik_hat;
vec3_scale(dfc*gijk*ex_delr,rik_hat,drk);
vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk);
vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk);
vec3_scale(prefactor,drk,drk);
}
diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp
index bdfc9ea15..291af44d8 100644
--- a/src/MEAM/pair_meam.cpp
+++ b/src/MEAM/pair_meam.cpp
@@ -1,946 +1,946 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) 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.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Greg Wagner (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_meam.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MAXLINE 1024
enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11,L12,B2};
int nkeywords = 21;
const char *keywords[] = {"Ec","alpha","rho0","delta","lattce",
"attrac","repuls","nn2","Cmin","Cmax","rc","delr",
"augt1","gsmooth_factor","re","ialloy",
"mixture_ref_t","erose_form","zbl",
"emb_lin_neg","bkgd_dyn"};
/* ---------------------------------------------------------------------- */
PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nmax = 0;
rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
maxneigh = 0;
scrfcn = dscrfcn = fcpair = NULL;
nelements = 0;
elements = NULL;
mass = NULL;
// set comm size needed by this Pair
comm_forward = 38;
comm_reverse = 30;
}
/* ----------------------------------------------------------------------
free all arrays
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairMEAM::~PairMEAM()
{
meam_cleanup_();
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
delete [] fmap;
}
}
/* ---------------------------------------------------------------------- */
void PairMEAM::compute(int eflag, int vflag)
{
int i,j,ii,n,inum_half,errorflag;
double evdwl;
int *ilist_half,*numneigh_half,**firstneigh_half;
int *numneigh_full,**firstneigh_full;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = eflag_global = vflag_global =
eflag_atom = vflag_atom = 0;
// grow local arrays if necessary
if (atom->nmax > nmax) {
memory->destroy(rho);
memory->destroy(rho0);
memory->destroy(rho1);
memory->destroy(rho2);
memory->destroy(rho3);
memory->destroy(frhop);
memory->destroy(gamma);
memory->destroy(dgamma1);
memory->destroy(dgamma2);
memory->destroy(dgamma3);
memory->destroy(arho2b);
memory->destroy(arho1);
memory->destroy(arho2);
memory->destroy(arho3);
memory->destroy(arho3b);
memory->destroy(t_ave);
memory->destroy(tsq_ave);
nmax = atom->nmax;
memory->create(rho,nmax,"pair:rho");
memory->create(rho0,nmax,"pair:rho0");
memory->create(rho1,nmax,"pair:rho1");
memory->create(rho2,nmax,"pair:rho2");
memory->create(rho3,nmax,"pair:rho3");
memory->create(frhop,nmax,"pair:frhop");
memory->create(gamma,nmax,"pair:gamma");
memory->create(dgamma1,nmax,"pair:dgamma1");
memory->create(dgamma2,nmax,"pair:dgamma2");
memory->create(dgamma3,nmax,"pair:dgamma3");
memory->create(arho2b,nmax,"pair:arho2b");
memory->create(arho1,nmax,3,"pair:arho1");
memory->create(arho2,nmax,6,"pair:arho2");
memory->create(arho3,nmax,10,"pair:arho3");
memory->create(arho3b,nmax,3,"pair:arho3b");
memory->create(t_ave,nmax,3,"pair:t_ave");
memory->create(tsq_ave,nmax,3,"pair:tsq_ave");
}
// neighbor list info
inum_half = listhalf->inum;
ilist_half = listhalf->ilist;
numneigh_half = listhalf->numneigh;
firstneigh_half = listhalf->firstneigh;
numneigh_full = listfull->numneigh;
firstneigh_full = listfull->firstneigh;
// strip neighbor lists of any special bond flags before using with MEAM
// necessary before doing neigh_f2c and neigh_c2f conversions each step
if (neighbor->ago == 0) {
neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full);
}
// check size of scrfcn based on half neighbor list
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
n = 0;
for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
if (n > maxneigh) {
memory->destroy(scrfcn);
memory->destroy(dscrfcn);
memory->destroy(fcpair);
maxneigh = n;
memory->create(scrfcn,maxneigh,"pair:scrfcn");
memory->create(dscrfcn,maxneigh,"pair:dscrfcn");
memory->create(fcpair,maxneigh,"pair:fcpair");
}
// zero out local arrays
for (i = 0; i < nall; i++) {
rho0[i] = 0.0;
arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
}
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int ntype = atom->ntypes;
// change neighbor list indices to Fortran indexing
neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full);
// 3 stages of MEAM calculation
// loop over my atoms followed by communication
int ifort;
int offset = 0;
errorflag = 0;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
ifort = i+1;
meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0],
&numneigh_half[i],firstneigh_half[i],
&numneigh_full[i],firstneigh_full[i],
&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
rho0,&arho1[0][0],&arho2[0][0],arho2b,
&arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],
&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
offset += numneigh_half[i];
}
comm->reverse_comm_pair(this);
meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
&eng_vdwl,eatom,&ntype,type,fmap,
&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
comm->forward_comm_pair(this);
offset = 0;
// vptr is first value in vatom if it will be used by meam_force()
// else vatom may not exist, so pass dummy ptr
double *vptr;
if (vflag_atom) vptr = &vatom[0][0];
else vptr = &cutmax;
for (ii = 0; ii < inum_half; ii++) {
i = ilist_half[ii];
ifort = i+1;
meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom,
&vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
&numneigh_half[i],firstneigh_half[i],
&numneigh_full[i],firstneigh_full[i],
&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
&t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->one(FLERR,str);
}
offset += numneigh_half[i];
}
// change neighbor list indices back to C indexing
neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half);
neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full);
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairMEAM::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
map = new int[n+1];
fmap = new int[n];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairMEAM::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairMEAM::coeff(int narg, char **arg)
{
int i,j,m,n;
if (!allocated) allocate();
if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read MEAM element names between 2 filenames
// nelements = # of MEAM elements
// elements = list of unique element names
if (nelements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
delete [] mass;
}
nelements = narg - 4 - atom->ntypes;
if (nelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
elements = new char*[nelements];
mass = new double[nelements];
for (i = 0; i < nelements; i++) {
n = strlen(arg[i+3]) + 1;
elements[i] = new char[n];
strcpy(elements[i],arg[i+3]);
}
// read MEAM library and parameter files
// pass all parameters to MEAM package
// tell MEAM package that setup is done
read_files(arg[2],arg[2+nelements+1]);
meam_setup_done_(&cutmax);
// read args that map atom types to MEAM elements
// map[i] = which element the Ith atom type is, -1 if not mapped
for (i = 4 + nelements; i < narg; i++) {
m = i - (4+nelements) + 1;
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
if (j < nelements) map[m] = j;
else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
else error->all(FLERR,"Incorrect args for pair coefficients");
}
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
// set mass for i,i in atom class
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
if (i == j) atom->set_mass(i,mass[map[i]]);
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairMEAM::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style MEAM requires newton pair on");
// need full and half neighbor list
int irequest_full = neighbor->request(this);
neighbor->requests[irequest_full]->id = 1;
neighbor->requests[irequest_full]->half = 0;
neighbor->requests[irequest_full]->full = 1;
int irequest_half = neighbor->request(this);
neighbor->requests[irequest_half]->id = 2;
neighbor->requests[irequest_half]->half = 0;
neighbor->requests[irequest_half]->half_from_full = 1;
neighbor->requests[irequest_half]->otherlist = irequest_full;
// setup Fortran-style mapping array needed by MEAM package
// fmap is indexed from 1:ntypes by Fortran and stores a Fortran index
// if type I is not a MEAM atom, fmap stores a 0
for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1;
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
half or full
------------------------------------------------------------------------- */
void PairMEAM::init_list(int id, NeighList *ptr)
{
if (id == 1) listfull = ptr;
else if (id == 2) listhalf = ptr;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairMEAM::init_one(int i, int j)
{
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairMEAM::read_files(char *globalfile, char *userfile)
{
// open global meamf file on proc 0
FILE *fp;
if (comm->me == 0) {
- fp = fopen(globalfile,"r");
+ fp = open_potential(globalfile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open MEAM potential file %s",globalfile);
error->one(FLERR,str);
}
}
// allocate parameter arrays
int params_per_line = 19;
int *lat = new int[nelements];
int *ielement = new int[nelements];
int *ibar = new int[nelements];
double *z = new double[nelements];
double *atwt = new double[nelements];
double *alpha = new double[nelements];
double *b0 = new double[nelements];
double *b1 = new double[nelements];
double *b2 = new double[nelements];
double *b3 = new double[nelements];
double *alat = new double[nelements];
double *esub = new double[nelements];
double *asub = new double[nelements];
double *t0 = new double[nelements];
double *t1 = new double[nelements];
double *t2 = new double[nelements];
double *t3 = new double[nelements];
double *rozero = new double[nelements];
bool *found = new bool[nelements];
for (int i = 0; i < nelements; i++) found[i] = false;
// read each set of params from global MEAM file
// one set of params can span multiple lines
// store params if element name is in element list
// if element name appears multiple times, only store 1st entry
int i,n,nwords;
char **words = new char*[params_per_line+1];
char line[MAXLINE],*ptr;
int eof = 0;
int nset = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if (ptr = strchr(line,'#')) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if (ptr = strchr(line,'#')) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in MEAM potential file");
// words = ptrs to all words in line
// strip single and double quotes from words
nwords = 0;
words[nwords++] = strtok(line,"' \t\n\r\f");
while (words[nwords++] = strtok(NULL,"' \t\n\r\f")) continue;
// skip if element name isn't in element list
for (i = 0; i < nelements; i++)
if (strcmp(words[0],elements[i]) == 0) break;
if (i == nelements) continue;
// skip if element already appeared
if (found[i] == true) continue;
found[i] = true;
// map lat string to an integer
if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 1");
// store parameters
z[i] = atof(words[2]);
ielement[i] = atoi(words[3]);
atwt[i] = atof(words[4]);
alpha[i] = atof(words[5]);
b0[i] = atof(words[6]);
b1[i] = atof(words[7]);
b2[i] = atof(words[8]);
b3[i] = atof(words[9]);
alat[i] = atof(words[10]);
esub[i] = atof(words[11]);
asub[i] = atof(words[12]);
t0[i] = atof(words[13]);
t1[i] = atof(words[14]);
t2[i] = atof(words[15]);
t3[i] = atof(words[16]);
rozero[i] = atof(words[17]);
ibar[i] = atoi(words[18]);
nset++;
}
// error if didn't find all elements in file
if (nset != nelements)
error->all(FLERR,"Did not find all elements in MEAM library file");
// pass element parameters to MEAM package
meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
// set element masses
for (i = 0; i < nelements; i++) mass[i] = atwt[i];
// clean-up memory
delete [] words;
delete [] lat;
delete [] ielement;
delete [] ibar;
delete [] z;
delete [] atwt;
delete [] alpha;
delete [] b0;
delete [] b1;
delete [] b2;
delete [] b3;
delete [] alat;
delete [] esub;
delete [] asub;
delete [] t0;
delete [] t1;
delete [] t2;
delete [] t3;
delete [] rozero;
delete [] found;
// done if user param file is NULL
if (strcmp(userfile,"NULL") == 0) return;
// open user param file on proc 0
if (comm->me == 0) {
- fp = fopen(userfile,"r");
+ fp = open_potential(userfile);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open MEAM potential file %s",userfile);
error->one(FLERR,str);
}
}
// read settings
// pass them one at a time to MEAM package
// match strings to list of corresponding ints
int which;
double value;
int nindex,index[3];
int maxparams = 6;
char **params = new char*[maxparams];
int nparams;
eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if (ptr = strchr(line,'#')) *ptr = '\0';
nparams = atom->count_words(line);
if (nparams == 0) continue;
// words = ptrs to all words in line
nparams = 0;
params[nparams++] = strtok(line,"=(), '\t\n\r\f");
while (nparams < maxparams &&
(params[nparams++] = strtok(NULL,"=(), '\t\n\r\f")))
continue;
nparams--;
for (which = 0; which < nkeywords; which++)
if (strcmp(params[0],keywords[which]) == 0) break;
if (which == nkeywords) {
char str[128];
sprintf(str,"Keyword %s in MEAM parameter file not recognized",
params[0]);
error->all(FLERR,str);
}
nindex = nparams - 2;
for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]);
// map lattce_meam value to an integer
if (which == 4) {
if (strcmp(params[nparams-1],"fcc") == 0) value = FCC;
else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC;
else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP;
else if (strcmp(params[nparams-1],"dim") == 0) value = DIM;
else if (strcmp(params[nparams-1],"dia") == 0) value = DIAMOND;
else if (strcmp(params[nparams-1],"b1") == 0) value = B1;
else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
else if (strcmp(params[nparams-1],"b2") == 0) value = B2;
else error->all(FLERR,"Unrecognized lattice type in MEAM file 2");
}
else value = atof(params[nparams-1]);
// pass single setting to MEAM package
int errorflag = 0;
meam_setup_param_(&which,&value,&nindex,index,&errorflag);
if (errorflag) {
char str[128];
sprintf(str,"MEAM library error %d",errorflag);
error->all(FLERR,str);
}
}
delete [] params;
}
/* ---------------------------------------------------------------------- */
int PairMEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
{
int i,j,k,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = rho0[j];
buf[m++] = rho1[j];
buf[m++] = rho2[j];
buf[m++] = rho3[j];
buf[m++] = frhop[j];
buf[m++] = gamma[j];
buf[m++] = dgamma1[j];
buf[m++] = dgamma2[j];
buf[m++] = dgamma3[j];
buf[m++] = arho2b[j];
buf[m++] = arho1[j][0];
buf[m++] = arho1[j][1];
buf[m++] = arho1[j][2];
buf[m++] = arho2[j][0];
buf[m++] = arho2[j][1];
buf[m++] = arho2[j][2];
buf[m++] = arho2[j][3];
buf[m++] = arho2[j][4];
buf[m++] = arho2[j][5];
for (k = 0; k < 10; k++) buf[m++] = arho3[j][k];
buf[m++] = arho3b[j][0];
buf[m++] = arho3b[j][1];
buf[m++] = arho3b[j][2];
buf[m++] = t_ave[j][0];
buf[m++] = t_ave[j][1];
buf[m++] = t_ave[j][2];
buf[m++] = tsq_ave[j][0];
buf[m++] = tsq_ave[j][1];
buf[m++] = tsq_ave[j][2];
}
return comm_forward;
}
/* ---------------------------------------------------------------------- */
void PairMEAM::unpack_comm(int n, int first, double *buf)
{
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
rho0[i] = buf[m++];
rho1[i] = buf[m++];
rho2[i] = buf[m++];
rho3[i] = buf[m++];
frhop[i] = buf[m++];
gamma[i] = buf[m++];
dgamma1[i] = buf[m++];
dgamma2[i] = buf[m++];
dgamma3[i] = buf[m++];
arho2b[i] = buf[m++];
arho1[i][0] = buf[m++];
arho1[i][1] = buf[m++];
arho1[i][2] = buf[m++];
arho2[i][0] = buf[m++];
arho2[i][1] = buf[m++];
arho2[i][2] = buf[m++];
arho2[i][3] = buf[m++];
arho2[i][4] = buf[m++];
arho2[i][5] = buf[m++];
for (k = 0; k < 10; k++) arho3[i][k] = buf[m++];
arho3b[i][0] = buf[m++];
arho3b[i][1] = buf[m++];
arho3b[i][2] = buf[m++];
t_ave[i][0] = buf[m++];
t_ave[i][1] = buf[m++];
t_ave[i][2] = buf[m++];
tsq_ave[i][0] = buf[m++];
tsq_ave[i][1] = buf[m++];
tsq_ave[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
{
int i,k,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
buf[m++] = rho0[i];
buf[m++] = arho2b[i];
buf[m++] = arho1[i][0];
buf[m++] = arho1[i][1];
buf[m++] = arho1[i][2];
buf[m++] = arho2[i][0];
buf[m++] = arho2[i][1];
buf[m++] = arho2[i][2];
buf[m++] = arho2[i][3];
buf[m++] = arho2[i][4];
buf[m++] = arho2[i][5];
for (k = 0; k < 10; k++) buf[m++] = arho3[i][k];
buf[m++] = arho3b[i][0];
buf[m++] = arho3b[i][1];
buf[m++] = arho3b[i][2];
buf[m++] = t_ave[i][0];
buf[m++] = t_ave[i][1];
buf[m++] = t_ave[i][2];
buf[m++] = tsq_ave[i][0];
buf[m++] = tsq_ave[i][1];
buf[m++] = tsq_ave[i][2];
}
return comm_reverse;
}
/* ---------------------------------------------------------------------- */
void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,k,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
rho0[j] += buf[m++];
arho2b[j] += buf[m++];
arho1[j][0] += buf[m++];
arho1[j][1] += buf[m++];
arho1[j][2] += buf[m++];
arho2[j][0] += buf[m++];
arho2[j][1] += buf[m++];
arho2[j][2] += buf[m++];
arho2[j][3] += buf[m++];
arho2[j][4] += buf[m++];
arho2[j][5] += buf[m++];
for (k = 0; k < 10; k++) arho3[j][k] += buf[m++];
arho3b[j][0] += buf[m++];
arho3b[j][1] += buf[m++];
arho3b[j][2] += buf[m++];
t_ave[j][0] += buf[m++];
t_ave[j][1] += buf[m++];
t_ave[j][2] += buf[m++];
tsq_ave[j][0] += buf[m++];
tsq_ave[j][1] += buf[m++];
tsq_ave[j][2] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double PairMEAM::memory_usage()
{
double bytes = 11 * nmax * sizeof(double);
bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double);
bytes += 3 * maxneigh * sizeof(double);
return bytes;
}
/* ----------------------------------------------------------------------
strip special bond flags from neighbor list entries
are not used with MEAM
need to do here so Fortran lib doesn't see them
done once per reneighbor so that neigh_f2c and neigh_c2f don't see them
------------------------------------------------------------------------- */
void PairMEAM::neigh_strip(int inum, int *ilist,
int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
}
}
/* ----------------------------------------------------------------------
toggle neighbor list indices between zero- and one-based values
needed for access by MEAM Fortran library
------------------------------------------------------------------------- */
void PairMEAM::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j]--;
}
}
void PairMEAM::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh)
{
int i,j,ii,jnum;
int *jlist;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i];
jnum = numneigh[i];
for (j = 0; j < jnum; j++) jlist[j]++;
}
}
diff --git a/src/USER-OMP/pair_tersoff_mod_omp.cpp b/src/USER-OMP/pair_tersoff_mod_omp.cpp
new file mode 100644
index 000000000..9a3f6f458
--- /dev/null
+++ b/src/USER-OMP/pair_tersoff_mod_omp.cpp
@@ -0,0 +1,252 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_tersoff_mod_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+
+#include "suffix.h"
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+PairTersoffMODOMP::PairTersoffMODOMP(LAMMPS *lmp) :
+ PairTersoffMOD(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ suffix_flag |= Suffix::OMP;
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairTersoffMODOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = vflag_atom = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ thr->timer(Timer::START);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (evflag) {
+ if (eflag) {
+ if (vflag_atom) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (vflag_atom) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else eval<0,0,0>(ifrom, ito, thr);
+
+ thr->timer(Timer::PAIR);
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
+void PairTersoffMODOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ int i,j,k,ii,jj,kk,jnum;
+ int itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+ double rsq,rsq1,rsq2;
+ double delr1[3],delr2[3],fi[3],fj[3],fk[3];
+ double zeta_ij,prefactor;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+
+ evdwl = 0.0;
+
+ const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
+ dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
+ const int * _noalias const tag = atom->tag;
+ const int * _noalias const type = atom->type;
+ const int nlocal = atom->nlocal;
+
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ double fxtmp,fytmp,fztmp;
+
+ // loop over full neighbor list of my atoms
+
+ for (ii = iifrom; ii < iito; ++ii) {
+
+ i = ilist[ii];
+ itag = tag[i];
+ itype = map[type[i]];
+ xtmp = x[i].x;
+ ytmp = x[i].y;
+ ztmp = x[i].z;
+ fxtmp = fytmp = fztmp = 0.0;
+
+ // two-body interactions, skip half of them
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ jtag = tag[j];
+
+ if (itag > jtag) {
+ if ((itag+jtag) % 2 == 0) continue;
+ } else if (itag < jtag) {
+ if ((itag+jtag) % 2 == 1) continue;
+ } else {
+ if (x[j].z < ztmp) continue;
+ if (x[j].z == ztmp && x[j].y < ytmp) continue;
+ if (x[j].z == ztmp && x[j].y == ytmp && x[j].x < xtmp) continue;
+ }
+
+ jtype = map[type[j]];
+
+ delx = xtmp - x[j].x;
+ dely = ytmp - x[j].y;
+ delz = ztmp - x[j].z;
+ rsq = delx*delx + dely*dely + delz*delz;
+
+ iparam_ij = elem2param[itype][jtype][jtype];
+ if (rsq > params[iparam_ij].cutsq) continue;
+
+ repulsive(¶ms[iparam_ij],rsq,fpair,EFLAG,evdwl);
+
+ fxtmp += delx*fpair;
+ fytmp += dely*fpair;
+ fztmp += delz*fpair;
+ f[j].x -= delx*fpair;
+ f[j].y -= dely*fpair;
+ f[j].z -= delz*fpair;
+
+ if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
+ evdwl,0.0,fpair,delx,dely,delz,thr);
+ }
+
+ // three-body interactions
+ // skip immediately if I-J is not within cutoff
+ double fjxtmp,fjytmp,fjztmp;
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+ jtype = map[type[j]];
+ iparam_ij = elem2param[itype][jtype][jtype];
+
+ delr1[0] = x[j].x - xtmp;
+ delr1[1] = x[j].y - ytmp;
+ delr1[2] = x[j].z - ztmp;
+ rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
+ if (rsq1 > params[iparam_ij].cutsq) continue;
+
+ // accumulate bondorder zeta for each i-j interaction via loop over k
+
+ fjxtmp = fjytmp = fjztmp = 0.0;
+ zeta_ij = 0.0;
+
+ for (kk = 0; kk < jnum; kk++) {
+ if (jj == kk) continue;
+ k = jlist[kk];
+ k &= NEIGHMASK;
+ ktype = map[type[k]];
+ iparam_ijk = elem2param[itype][jtype][ktype];
+
+ delr2[0] = x[k].x - xtmp;
+ delr2[1] = x[k].y - ytmp;
+ delr2[2] = x[k].z - ztmp;
+ rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
+ if (rsq2 > params[iparam_ijk].cutsq) continue;
+
+ zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
+ }
+
+ // pairwise force due to zeta
+
+ force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fpair,prefactor,EFLAG,evdwl);
+
+ fxtmp += delr1[0]*fpair;
+ fytmp += delr1[1]*fpair;
+ fztmp += delr1[2]*fpair;
+ fjxtmp -= delr1[0]*fpair;
+ fjytmp -= delr1[1]*fpair;
+ fjztmp -= delr1[2]*fpair;
+
+ if (EVFLAG) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,evdwl,0.0,
+ -fpair,-delr1[0],-delr1[1],-delr1[2],thr);
+
+ // attractive term via loop over k
+
+ for (kk = 0; kk < jnum; kk++) {
+ if (jj == kk) continue;
+ k = jlist[kk];
+ k &= NEIGHMASK;
+ ktype = map[type[k]];
+ iparam_ijk = elem2param[itype][jtype][ktype];
+
+ delr2[0] = x[k].x - xtmp;
+ delr2[1] = x[k].y - ytmp;
+ delr2[2] = x[k].z - ztmp;
+ rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
+ if (rsq2 > params[iparam_ijk].cutsq) continue;
+
+ attractive(¶ms[iparam_ijk],prefactor,
+ rsq1,rsq2,delr1,delr2,fi,fj,fk);
+
+ fxtmp += fi[0];
+ fytmp += fi[1];
+ fztmp += fi[2];
+ fjxtmp += fj[0];
+ fjytmp += fj[1];
+ fjztmp += fj[2];
+ f[k].x += fk[0];
+ f[k].y += fk[1];
+ f[k].z += fk[2];
+
+ if (VFLAG_ATOM) v_tally3_thr(i,j,k,fj,fk,delr1,delr2,thr);
+ }
+ f[j].x += fjxtmp;
+ f[j].y += fjytmp;
+ f[j].z += fjztmp;
+ }
+ f[i].x += fxtmp;
+ f[i].y += fytmp;
+ f[i].z += fztmp;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairTersoffMODOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairTersoffMOD::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_tersoff_mod_omp.h b/src/USER-OMP/pair_tersoff_mod_omp.h
new file mode 100644
index 000000000..ed0109449
--- /dev/null
+++ b/src/USER-OMP/pair_tersoff_mod_omp.h
@@ -0,0 +1,43 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(tersoff/mod/omp,PairTersoffMODOMP)
+
+#else
+
+#ifndef LMP_PAIR_TERSOFF_MOD_OMP_H
+#define LMP_PAIR_TERSOFF_MOD_OMP_H
+
+#include "pair_tersoff_mod.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairTersoffMODOMP : public PairTersoffMOD, public ThrOMP {
+
+ public:
+ PairTersoffMODOMP(class LAMMPS *);
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template <int EVFLAG, int EFLAG, int VFLAG_ATOM>
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/pair_zbl_omp.cpp b/src/USER-OMP/pair_zbl_omp.cpp
new file mode 100644
index 000000000..ad028b09a
--- /dev/null
+++ b/src/USER-OMP/pair_zbl_omp.cpp
@@ -0,0 +1,172 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ This software is distributed under the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "pair_zbl_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+
+#include "suffix.h"
+using namespace LAMMPS_NS;
+using namespace PairZBLConstants;
+
+/* ---------------------------------------------------------------------- */
+
+PairZBLOMP::PairZBLOMP(LAMMPS *lmp) :
+ PairZBL(lmp), ThrOMP(lmp, THR_PAIR)
+{
+ suffix_flag |= Suffix::OMP;
+ respa_enable = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairZBLOMP::compute(int eflag, int vflag)
+{
+ if (eflag || vflag) {
+ ev_setup(eflag,vflag);
+ } else evflag = vflag_fdotr = 0;
+
+ const int nall = atom->nlocal + atom->nghost;
+ const int nthreads = comm->nthreads;
+ const int inum = list->inum;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+ {
+ int ifrom, ito, tid;
+
+ loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+ ThrData *thr = fix->get_thr(tid);
+ thr->timer(Timer::START);
+ ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+ if (evflag) {
+ if (eflag) {
+ if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr);
+ else eval<1,1,0>(ifrom, ito, thr);
+ } else {
+ if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr);
+ else eval<1,0,0>(ifrom, ito, thr);
+ }
+ } else {
+ if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr);
+ else eval<0,0,0>(ifrom, ito, thr);
+ }
+ thr->timer(Timer::PAIR);
+ reduce_thr(this, eflag, vflag, thr);
+ } // end of omp parallel region
+}
+
+template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
+void PairZBLOMP::eval(int iifrom, int iito, ThrData * const thr)
+{
+ const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
+ dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
+ const int * _noalias const type = atom->type;
+ const int * _noalias const ilist = list->ilist;
+ const int * _noalias const numneigh = list->numneigh;
+ const int * const * const firstneigh = list->firstneigh;
+
+ double xtmp,ytmp,ztmp,delx,dely,delz,fxtmp,fytmp,fztmp;
+ double rsq,t,fswitch,eswitch,evdwl,fpair;
+
+ const int nlocal = atom->nlocal;
+ int j,jj,jnum,jtype;
+
+ evdwl = 0.0;
+
+ // loop over neighbors of my atoms
+
+ for (int ii = iifrom; ii < iito; ++ii) {
+ const int i = ilist[ii];
+ const int itype = type[i];
+ const int * _noalias const jlist = firstneigh[i];
+ const double * _noalias const sw1i = sw1[itype];
+ const double * _noalias const sw2i = sw2[itype];
+ const double * _noalias const sw3i = sw3[itype];
+ const double * _noalias const sw4i = sw4[itype];
+ const double * _noalias const sw5i = sw5[itype];
+
+ xtmp = x[i].x;
+ ytmp = x[i].y;
+ ztmp = x[i].z;
+ jnum = numneigh[i];
+ fxtmp=fytmp=fztmp=0.0;
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ j &= NEIGHMASK;
+
+ delx = xtmp - x[j].x;
+ dely = ytmp - x[j].y;
+ delz = ztmp - x[j].z;
+ rsq = delx*delx + dely*dely + delz*delz;
+ jtype = type[j];
+
+ if (rsq < cut_globalsq) {
+ const double r = sqrt(rsq);
+ fpair = dzbldr(r, itype, jtype);
+
+ if (r > cut_inner) {
+ t = r - cut_inner;
+ fswitch = t*t *
+ (sw1i[jtype] + sw2i[jtype]*t);
+ fpair += fswitch;
+ }
+
+ fpair *= -1.0/r;
+ fxtmp += delx*fpair;
+ fytmp += dely*fpair;
+ fztmp += delz*fpair;
+
+ if (NEWTON_PAIR || j < nlocal) {
+ f[j].x -= delx*fpair;
+ f[j].y -= dely*fpair;
+ f[j].z -= delz*fpair;
+ }
+
+ if (EFLAG) {
+ evdwl = e_zbl(r, itype, jtype);
+ evdwl += sw5i[jtype];
+ if (r > cut_inner) {
+ eswitch = t*t*t *
+ (sw3i[jtype] + sw4i[jtype]*t);
+ evdwl += eswitch;
+ }
+ }
+
+ if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR,
+ evdwl,0.0,fpair,delx,dely,delz,thr);
+ }
+ }
+ f[i].x += fxtmp;
+ f[i].y += fytmp;
+ f[i].z += fztmp;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairZBLOMP::memory_usage()
+{
+ double bytes = memory_usage_thr();
+ bytes += PairZBL::memory_usage();
+
+ return bytes;
+}
diff --git a/src/USER-OMP/pair_zbl_omp.h b/src/USER-OMP/pair_zbl_omp.h
new file mode 100644
index 000000000..a75d9dba5
--- /dev/null
+++ b/src/USER-OMP/pair_zbl_omp.h
@@ -0,0 +1,48 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) 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.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(zbl/omp,PairZBLOMP)
+
+#else
+
+#ifndef LMP_PAIR_ZBL_OMP_H
+#define LMP_PAIR_ZBL_OMP_H
+
+#include "pair_zbl.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class PairZBLOMP : public PairZBL, public ThrOMP {
+
+ public:
+ PairZBLOMP(class LAMMPS *);
+
+ virtual void compute(int, int);
+ virtual double memory_usage();
+
+ private:
+ template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
+ void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif
diff --git a/src/USER-OMP/thr_data.cpp b/src/USER-OMP/thr_data.cpp
index 3b3b87dbc..407c3fdea 100644
--- a/src/USER-OMP/thr_data.cpp
+++ b/src/USER-OMP/thr_data.cpp
@@ -1,310 +1,363 @@
/* -------------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) 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.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
per-thread data management for LAMMPS
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "thr_data.h"
#include <string.h>
#include <stdio.h>
#include "memory.h"
#include "timer.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ThrData::ThrData(int tid, Timer *t)
: _f(0),_torque(0),_erforce(0),_de(0),_drho(0),_mu(0),_lambda(0),_rhoB(0),
_D_values(0),_rho(0),_fp(0),_rho1d(0),_drho1d(0),_tid(tid), _timer(t)
{
_timer_active = 0;
}
/* ---------------------------------------------------------------------- */
void ThrData::check_tid(int tid)
{
if (tid != _tid)
fprintf(stderr,"WARNING: external and internal tid mismatch %d != %d\n",tid,_tid);
}
/* ---------------------------------------------------------------------- */
void ThrData::_stamp(enum Timer::ttype flag)
{
// do nothing until it gets set to 0 in ::setup()
if (_timer_active < 0) return;
if (flag == Timer::START) {
_timer_active = 1;
}
if (_timer_active) _timer->stamp(flag);
}
/* ---------------------------------------------------------------------- */
double ThrData::get_time(enum Timer::ttype flag)
{
if (_timer)
return _timer->get_wall(flag);
else
return 0.0;
}
/* ---------------------------------------------------------------------- */
void ThrData::init_force(int nall, double **f, double **torque,
double *erforce, double *de, double *drho)
{
eng_vdwl=eng_coul=eng_bond=eng_angle=eng_dihed=eng_imprp=eng_kspce=0.0;
memset(virial_pair,0,6*sizeof(double));
memset(virial_bond,0,6*sizeof(double));
memset(virial_angle,0,6*sizeof(double));
memset(virial_dihed,0,6*sizeof(double));
memset(virial_imprp,0,6*sizeof(double));
memset(virial_kspce,0,6*sizeof(double));
eatom_pair=eatom_bond=eatom_angle=eatom_dihed=eatom_imprp=eatom_kspce=NULL;
vatom_pair=vatom_bond=vatom_angle=vatom_dihed=vatom_imprp=vatom_kspce=NULL;
_f = f + _tid*nall;
if (nall > 0)
memset(&(_f[0][0]),0,nall*3*sizeof(double));
if (torque) {
_torque = torque + _tid*nall;
if (nall > 0)
memset(&(_torque[0][0]),0,nall*3*sizeof(double));
} else _torque = NULL;
if (erforce) {
_erforce = erforce + _tid*nall;
if (nall > 0)
memset(&(_erforce[0]),0,nall*sizeof(double));
} else _erforce = NULL;
if (de) {
_de = de + _tid*nall;
if (nall > 0)
memset(&(_de[0]),0,nall*sizeof(double));
} else _de = NULL;
if (drho) {
_drho = drho + _tid*nall;
if (nall > 0)
memset(&(_drho[0]),0,nall*sizeof(double));
} else _drho = NULL;
}
/* ----------------------------------------------------------------------
set up and clear out locally managed per atom arrays
------------------------------------------------------------------------- */
void ThrData::init_eam(int nall, double *rho)
{
_rho = rho + _tid*nall;
if (nall > 0)
memset(_rho, 0, nall*sizeof(double));
}
/* ---------------------------------------------------------------------- */
void ThrData::init_adp(int nall, double *rho, double **mu, double **lambda)
{
init_eam(nall, rho);
_mu = mu + _tid*nall;
_lambda = lambda + _tid*nall;
if (nall > 0) {
memset(&(_mu[0][0]), 0, nall*3*sizeof(double));
memset(&(_lambda[0][0]), 0, nall*6*sizeof(double));
}
}
/* ---------------------------------------------------------------------- */
void ThrData::init_cdeam(int nall, double *rho, double *rhoB, double *D_values)
{
init_eam(nall, rho);
_rhoB = rhoB + _tid*nall;
_D_values = D_values + _tid*nall;
if (nall > 0) {
memset(_rhoB, 0, nall*sizeof(double));
memset(_D_values, 0, nall*sizeof(double));
}
}
/* ---------------------------------------------------------------------- */
void ThrData::init_eim(int nall, double *rho, double *fp)
{
init_eam(nall, rho);
_fp = fp + _tid*nall;
if (nall > 0)
memset(_fp,0,nall*sizeof(double));
}
/* ----------------------------------------------------------------------
if order > 0 : set up per thread storage for PPPM
if order < 0 : free per thread storage for PPPM
------------------------------------------------------------------------- */
#if defined(FFT_SINGLE)
typedef float FFT_SCALAR;
#else
typedef double FFT_SCALAR;
#endif
void ThrData::init_pppm(int order, Memory *memory)
{
FFT_SCALAR **rho1d, **drho1d;
if (order > 0) {
memory->create2d_offset(rho1d,3,-order/2,order/2,"thr_data:rho1d");
memory->create2d_offset(drho1d,3,-order/2,order/2,"thr_data:drho1d");
_rho1d = static_cast<void *>(rho1d);
_drho1d = static_cast<void *>(drho1d);
} else {
order = -order;
rho1d = static_cast<FFT_SCALAR **>(_rho1d);
drho1d = static_cast<FFT_SCALAR **>(_drho1d);
memory->destroy2d_offset(rho1d,-order/2);
memory->destroy2d_offset(drho1d,-order/2);
}
}
/* ----------------------------------------------------------------------
if order > 0 : set up per thread storage for PPPM
if order < 0 : free per thread storage for PPPM
------------------------------------------------------------------------- */
#if defined(FFT_SINGLE)
typedef float FFT_SCALAR;
#else
typedef double FFT_SCALAR;
#endif
void ThrData::init_pppm_disp(int order_6, Memory *memory)
{
FFT_SCALAR **rho1d_6, **drho1d_6;
if (order_6 > 0) {
memory->create2d_offset(rho1d_6,3,-order_6/2,order_6/2,"thr_data:rho1d_6");
memory->create2d_offset(drho1d_6,3,-order_6/2,order_6/2,"thr_data:drho1d_6");
_rho1d_6 = static_cast<void *>(rho1d_6);
_drho1d_6 = static_cast<void *>(drho1d_6);
} else {
order_6 = -order_6;
rho1d_6 = static_cast<FFT_SCALAR **>(_rho1d_6);
drho1d_6 = static_cast<FFT_SCALAR **>(_drho1d_6);
memory->destroy2d_offset(rho1d_6,-order_6/2);
memory->destroy2d_offset(drho1d_6,-order_6/2);
}
}
/* ----------------------------------------------------------------------
compute global pair virial via summing F dot r over own & ghost atoms
at this point, only pairwise forces have been accumulated in atom->f
------------------------------------------------------------------------- */
void ThrData::virial_fdotr_compute(double **x, int nlocal, int nghost, int nfirst)
{
// sum over force on all particles including ghosts
if (nfirst < 0) {
int nall = nlocal + nghost;
for (int i = 0; i < nall; i++) {
virial_pair[0] += _f[i][0]*x[i][0];
virial_pair[1] += _f[i][1]*x[i][1];
virial_pair[2] += _f[i][2]*x[i][2];
virial_pair[3] += _f[i][1]*x[i][0];
virial_pair[4] += _f[i][2]*x[i][0];
virial_pair[5] += _f[i][2]*x[i][1];
}
// neighbor includegroup flag is set
// sum over force on initial nfirst particles and ghosts
} else {
int nall = nfirst;
for (int i = 0; i < nall; i++) {
virial_pair[0] += _f[i][0]*x[i][0];
virial_pair[1] += _f[i][1]*x[i][1];
virial_pair[2] += _f[i][2]*x[i][2];
virial_pair[3] += _f[i][1]*x[i][0];
virial_pair[4] += _f[i][2]*x[i][0];
virial_pair[5] += _f[i][2]*x[i][1];
}
nall = nlocal + nghost;
for (int i = nlocal; i < nall; i++) {
virial_pair[0] += _f[i][0]*x[i][0];
virial_pair[1] += _f[i][1]*x[i][1];
virial_pair[2] += _f[i][2]*x[i][2];
virial_pair[3] += _f[i][1]*x[i][0];
virial_pair[4] += _f[i][2]*x[i][0];
virial_pair[5] += _f[i][2]*x[i][1];
}
}
}
/* ---------------------------------------------------------------------- */
double ThrData::memory_usage()
{
double bytes = (7 + 6*6) * sizeof(double);
bytes += 2 * sizeof(double*);
bytes += 4 * sizeof(int);
return bytes;
}
/* additional helper functions */
// reduce per thread data into the first part of the data
// array that is used for the non-threaded parts and reset
// the temporary storage to 0.0. this routine depends on
// multi-dimensional arrays like force stored in this order
// x1,y1,z1,x2,y2,z2,...
// we need to post a barrier to wait until all threads are done
// with writing to the array .
void LAMMPS_NS::data_reduce_thr(double *dall, int nall, int nthreads, int ndim, int tid)
{
#if defined(_OPENMP)
// NOOP in single-threaded execution.
if (nthreads == 1) return;
#pragma omp barrier
{
const int nvals = ndim*nall;
const int idelta = nvals/nthreads + 1;
const int ifrom = tid*idelta;
const int ito = ((ifrom + idelta) > nvals) ? nvals : (ifrom + idelta);
- for (int m = ifrom; m < ito; ++m) {
- for (int n = 1; n < nthreads; ++n) {
- dall[m] += dall[n*nvals + m];
- dall[n*nvals + m] = 0.0;
+#if defined(USER_OMP_NO_UNROLL)
+ if (ifrom < nvals) {
+ int m = 0;
+
+ for (m = ifrom; m < ito; ++m) {
+ for (int n = 1; n < nthreads; ++n) {
+ dall[m] += dall[n*nvals + m];
+ dall[n*nvals + m] = 0.0;
+ }
}
}
+#else
+ // this if protects against having more threads than atoms
+ if (ifrom < nvals) {
+ int m = 0;
+
+ // for architectures that have L1 D-cache line sizes of 64 bytes
+ // (8 doubles) wide, explictly unroll this loop to compute 8
+ // contiguous values in the array at a time
+ // -- modify this code based on the size of the cache line
+ double t0, t1, t2, t3, t4, t5, t6, t7;
+ for (m = ifrom; m < (ito-7); m+=8) {
+ t0 = dall[m ];
+ t1 = dall[m+1];
+ t2 = dall[m+2];
+ t3 = dall[m+3];
+ t4 = dall[m+4];
+ t5 = dall[m+5];
+ t6 = dall[m+6];
+ t7 = dall[m+7];
+ for (int n = 1; n < nthreads; ++n) {
+ t0 += dall[n*nvals + m ];
+ t1 += dall[n*nvals + m+1];
+ t2 += dall[n*nvals + m+2];
+ t3 += dall[n*nvals + m+3];
+ t4 += dall[n*nvals + m+4];
+ t5 += dall[n*nvals + m+5];
+ t6 += dall[n*nvals + m+6];
+ t7 += dall[n*nvals + m+7];
+ dall[n*nvals + m ] = 0.0;
+ dall[n*nvals + m+1] = 0.0;
+ dall[n*nvals + m+2] = 0.0;
+ dall[n*nvals + m+3] = 0.0;
+ dall[n*nvals + m+4] = 0.0;
+ dall[n*nvals + m+5] = 0.0;
+ dall[n*nvals + m+6] = 0.0;
+ dall[n*nvals + m+7] = 0.0;
+ }
+ dall[m ] = t0;
+ dall[m+1] = t1;
+ dall[m+2] = t2;
+ dall[m+3] = t3;
+ dall[m+4] = t4;
+ dall[m+5] = t5;
+ dall[m+6] = t6;
+ dall[m+7] = t7;
+ }
+ }
+#endif
}
#else
// NOOP in non-threaded execution.
return;
#endif
}
diff --git a/src/input.cpp b/src/input.cpp
index 17d32f7d7..1e40d1e4f 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -1,1578 +1,1591 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) 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.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "mpi.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "ctype.h"
#include "unistd.h"
#include "sys/stat.h"
#include "input.h"
#include "style_command.h"
#include "universe.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "group.h"
#include "domain.h"
#include "output.h"
#include "thermo.h"
#include "force.h"
#include "pair.h"
#include "min.h"
#include "modify.h"
#include "compute.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "update.h"
#include "neighbor.h"
#include "special.h"
+#include "timer.h"
#include "variable.h"
#include "accelerator_cuda.h"
#include "error.h"
#include "memory.h"
#ifdef _OPENMP
#include "omp.h"
#endif
using namespace LAMMPS_NS;
#define DELTALINE 256
#define DELTA 4
/* ---------------------------------------------------------------------- */
Input::Input(LAMMPS *lmp, int argc, char **argv) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
maxline = maxcopy = maxwork = 0;
line = copy = work = NULL;
narg = maxarg = 0;
arg = NULL;
echo_screen = 0;
echo_log = 1;
label_active = 0;
labelstr = NULL;
jump_skip = 0;
ifthenelse_flag = 0;
if (me == 0) {
nfile = maxfile = 1;
infiles = (FILE **) memory->smalloc(sizeof(FILE *),"input:infiles");
infiles[0] = infile;
} else infiles = NULL;
variable = new Variable(lmp);
// fill map with commands listed in style_command.h
command_map = new std::map<std::string,CommandCreator>();
#define COMMAND_CLASS
#define CommandStyle(key,Class) \
(*command_map)[#key] = &command_creator<Class>;
#include "style_command.h"
#undef CommandStyle
#undef COMMAND_CLASS
// process command-line args
// check for args "-var" and "-echo"
// caller has already checked that sufficient arguments exist
int iarg = 0;
while (iarg < argc) {
if (strcmp(argv[iarg],"-var") == 0 || strcmp(argv[iarg],"-v") == 0) {
int jarg = iarg+3;
while (jarg < argc && argv[jarg][0] != '-') jarg++;
variable->set(argv[iarg+1],jarg-iarg-2,&argv[iarg+2]);
iarg = jarg;
} else if (strcmp(argv[iarg],"-echo") == 0 ||
strcmp(argv[iarg],"-e") == 0) {
narg = 1;
char **tmp = arg; // trick echo() into using argv instead of arg
arg = &argv[iarg+1];
echo();
arg = tmp;
iarg += 2;
} else iarg++;
}
}
/* ---------------------------------------------------------------------- */
Input::~Input()
{
// don't free command and arg strings
// they just point to other allocated memory
memory->sfree(line);
memory->sfree(copy);
memory->sfree(work);
if (labelstr) delete [] labelstr;
memory->sfree(arg);
memory->sfree(infiles);
delete variable;
delete command_map;
}
/* ----------------------------------------------------------------------
process all input from infile
infile = stdin or file if command-line arg "-in" was used
------------------------------------------------------------------------- */
void Input::file()
{
int m,n;
while (1) {
// read a line from input script
// n = length of line including str terminator, 0 if end of file
// if line ends in continuation char '&', concatenate next line
if (me == 0) {
m = 0;
while (1) {
if (maxline-m < 2) reallocate(line,maxline,0);
if (fgets(&line[m],maxline-m,infile) == NULL) {
if (m) n = strlen(line) + 1;
else n = 0;
break;
}
m = strlen(line);
if (line[m-1] != '\n') continue;
m--;
while (m >= 0 && isspace(line[m])) m--;
if (m < 0 || line[m] != '&') {
line[m+1] = '\0';
n = m+2;
break;
}
}
}
// bcast the line
// if n = 0, end-of-file
// error if label_active is set, since label wasn't encountered
// if original input file, code is done
// else go back to previous input file
MPI_Bcast(&n,1,MPI_INT,0,world);
if (n == 0) {
if (label_active) error->all(FLERR,"Label wasn't found in input script");
if (me == 0) {
if (infile != stdin) fclose(infile);
nfile--;
}
MPI_Bcast(&nfile,1,MPI_INT,0,world);
if (nfile == 0) break;
if (me == 0) infile = infiles[nfile-1];
continue;
}
if (n > maxline) reallocate(line,maxline,n);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// echo the command unless scanning for label
if (me == 0 && label_active == 0) {
if (echo_screen && screen) fprintf(screen,"%s\n",line);
if (echo_log && logfile) fprintf(logfile,"%s\n",line);
}
// parse the line
// if no command, skip to next line in input script
parse();
if (command == NULL) continue;
// if scanning for label, skip command unless it's a label command
if (label_active && strcmp(command,"label") != 0) continue;
// execute the command
if (execute_command()) {
char *str = new char[maxline+32];
sprintf(str,"Unknown command: %s",line);
error->all(FLERR,str);
}
}
}
/* ----------------------------------------------------------------------
process all input from filename
called from library interface
------------------------------------------------------------------------- */
void Input::file(const char *filename)
{
// error if another nested file still open, should not be possible
// open new filename and set infile, infiles[0], nfile
// call to file() will close filename and decrement nfile
if (me == 0) {
if (nfile > 1)
error->one(FLERR,"Invalid use of library file() function");
infile = fopen(filename,"r");
if (infile == NULL) {
char str[128];
sprintf(str,"Cannot open input script %s",filename);
error->one(FLERR,str);
}
infiles[0] = infile;
nfile = 1;
}
file();
}
/* ----------------------------------------------------------------------
copy command in single to line, parse and execute it
return command name to caller
------------------------------------------------------------------------- */
char *Input::one(const char *single)
{
int n = strlen(single) + 1;
if (n > maxline) reallocate(line,maxline,n);
strcpy(line,single);
// echo the command unless scanning for label
if (me == 0 && label_active == 0) {
if (echo_screen && screen) fprintf(screen,"%s\n",line);
if (echo_log && logfile) fprintf(logfile,"%s\n",line);
}
// parse the line
// if no command, just return NULL
parse();
if (command == NULL) return NULL;
// if scanning for label, skip command unless it's a label command
if (label_active && strcmp(command,"label") != 0) return NULL;
// execute the command and return its name
if (execute_command()) {
char *str = new char[maxline+32];
sprintf(str,"Unknown command: %s",line);
error->all(FLERR,str);
}
return command;
}
/* ----------------------------------------------------------------------
parse copy of command line by inserting string terminators
strip comment = all chars from # on
replace all $ via variable substitution
command = first word
narg = # of args
arg[] = individual args
treat text between single/double quotes as one arg
------------------------------------------------------------------------- */
void Input::parse()
{
// duplicate line into copy string to break into words
int n = strlen(line) + 1;
if (n > maxcopy) reallocate(copy,maxcopy,n);
strcpy(copy,line);
// strip any # comment by replacing it with 0
// do not strip # inside single/double quotes
char quote = '\0';
char *ptr = copy;
while (*ptr) {
if (*ptr == '#' && !quote) {
*ptr = '\0';
break;
}
if (*ptr == quote) quote = '\0';
else if (*ptr == '"' || *ptr == '\'') quote = *ptr;
ptr++;
}
// perform $ variable substitution (print changes)
// except if searching for a label since earlier variable may not be defined
if (!label_active) substitute(copy,work,maxcopy,maxwork,1);
// command = 1st arg in copy string
char *next;
command = nextword(copy,&next);
if (command == NULL) return;
// point arg[] at each subsequent arg in copy string
// nextword() inserts string terminators into copy string to delimit args
// nextword() treats text between single/double quotes as one arg
narg = 0;
ptr = next;
while (ptr) {
if (narg == maxarg) {
maxarg += DELTA;
arg = (char **) memory->srealloc(arg,maxarg*sizeof(char *),"input:arg");
}
arg[narg] = nextword(ptr,&next);
if (!arg[narg]) break;
narg++;
ptr = next;
}
}
/* ----------------------------------------------------------------------
find next word in str
insert 0 at end of word
ignore leading whitespace
treat text between single/double quotes as one arg
matching quote must be followed by whitespace char if not end of string
strip quotes from returned word
return ptr to start of word
return next = ptr after word or NULL if word ended with 0
return NULL if no word in string
------------------------------------------------------------------------- */
char *Input::nextword(char *str, char **next)
{
char *start,*stop;
start = &str[strspn(str," \t\n\v\f\r")];
if (*start == '\0') return NULL;
if (*start == '"' || *start == '\'') {
stop = strchr(&start[1],*start);
if (!stop) error->all(FLERR,"Unbalanced quotes in input line");
if (stop[1] && !isspace(stop[1]))
error->all(FLERR,"Input line quote not followed by whitespace");
start++;
} else stop = &start[strcspn(start," \t\n\v\f\r")];
if (*stop == '\0') *next = NULL;
else *next = stop+1;
*stop = '\0';
return start;
}
/* ----------------------------------------------------------------------
substitute for $ variables in str using work str2 and return it
reallocate str/str2 to hold expanded version if necessary & reset max/max2
print updated string if flag is set and not searching for label
label_active will be 0 if called from external class
------------------------------------------------------------------------- */
void Input::substitute(char *&str, char *&str2, int &max, int &max2, int flag)
{
// use str2 as scratch space to expand str, then copy back to str
// reallocate str and str2 as necessary
// do not replace $ inside single/double quotes
// var = pts at variable name, ended by NULL
// if $ is followed by '{', trailing '}' becomes NULL
// else $x becomes x followed by NULL
// beyond = points to text following variable
int i,n,paren_count;
char immediate[256];
char *var,*value,*beyond;
char quote = '\0';
char *ptr = str;
n = strlen(str) + 1;
if (n > max2) reallocate(str2,max2,n);
*str2 = '\0';
char *ptr2 = str2;
while (*ptr) {
// variable substitution
if (*ptr == '$' && !quote) {
// value = ptr to expanded variable
// variable name between curly braces, e.g. ${a}
if (*(ptr+1) == '{') {
var = ptr+2;
i = 0;
while (var[i] != '\0' && var[i] != '}') i++;
if (var[i] == '\0') error->one(FLERR,"Invalid variable name");
var[i] = '\0';
beyond = ptr + strlen(var) + 3;
value = variable->retrieve(var);
// immediate variable between parenthesis, e.g. $(1/2)
} else if (*(ptr+1) == '(') {
var = ptr+2;
paren_count = 0;
i = 0;
while (var[i] != '\0' && !(var[i] == ')' && paren_count == 0)) {
switch (var[i]) {
case '(': paren_count++; break;
case ')': paren_count--; break;
default: ;
}
i++;
}
if (var[i] == '\0') error->one(FLERR,"Invalid immediate variable");
var[i] = '\0';
beyond = ptr + strlen(var) + 3;
sprintf(immediate,"%.20g",variable->compute_equal(var));
value = immediate;
// single character variable name, e.g. $a
} else {
var = ptr;
var[0] = var[1];
var[1] = '\0';
beyond = ptr + 2;
value = variable->retrieve(var);
}
if (value == NULL) error->one(FLERR,"Substitution for illegal variable");
// check if storage in str2 needs to be expanded
// re-initialize ptr and ptr2 to the point beyond the variable.
n = strlen(str2) + strlen(value) + strlen(beyond) + 1;
if (n > max2) reallocate(str2,max2,n);
strcat(str2,value);
ptr2 = str2 + strlen(str2);
ptr = beyond;
// output substitution progress if requested
if (flag && me == 0 && label_active == 0) {
if (echo_screen && screen) fprintf(screen,"%s%s\n",str2,beyond);
if (echo_log && logfile) fprintf(logfile,"%s%s\n",str2,beyond);
}
continue;
}
if (*ptr == quote) quote = '\0';
else if (*ptr == '"' || *ptr == '\'') quote = *ptr;
// copy current character into str2
*ptr2++ = *ptr++;
*ptr2 = '\0';
}
// set length of input str to length of work str2
// copy work string back to input str
if (max2 > max) reallocate(str,max,max2);
strcpy(str,str2);
}
/* ----------------------------------------------------------------------
rellocate a string
if n > 0: set max >= n in increments of DELTALINE
if n = 0: just increment max by DELTALINE
------------------------------------------------------------------------- */
void Input::reallocate(char *&str, int &max, int n)
{
if (n) {
while (n > max) max += DELTALINE;
} else max += DELTALINE;
str = (char *) memory->srealloc(str,max*sizeof(char),"input:str");
}
/* ----------------------------------------------------------------------
process a single parsed command
return 0 if successful, -1 if did not recognize command
------------------------------------------------------------------------- */
int Input::execute_command()
{
int flag = 1;
if (!strcmp(command,"clear")) clear();
else if (!strcmp(command,"echo")) echo();
else if (!strcmp(command,"if")) ifthenelse();
else if (!strcmp(command,"include")) include();
else if (!strcmp(command,"jump")) jump();
else if (!strcmp(command,"label")) label();
else if (!strcmp(command,"log")) log();
else if (!strcmp(command,"next")) next_command();
else if (!strcmp(command,"partition")) partition();
else if (!strcmp(command,"print")) print();
else if (!strcmp(command,"quit")) quit();
else if (!strcmp(command,"shell")) shell();
else if (!strcmp(command,"variable")) variable_command();
else if (!strcmp(command,"angle_coeff")) angle_coeff();
else if (!strcmp(command,"angle_style")) angle_style();
else if (!strcmp(command,"atom_modify")) atom_modify();
else if (!strcmp(command,"atom_style")) atom_style();
else if (!strcmp(command,"bond_coeff")) bond_coeff();
else if (!strcmp(command,"bond_style")) bond_style();
else if (!strcmp(command,"boundary")) boundary();
else if (!strcmp(command,"box")) box();
else if (!strcmp(command,"communicate")) communicate();
else if (!strcmp(command,"compute")) compute();
else if (!strcmp(command,"compute_modify")) compute_modify();
else if (!strcmp(command,"dielectric")) dielectric();
else if (!strcmp(command,"dihedral_coeff")) dihedral_coeff();
else if (!strcmp(command,"dihedral_style")) dihedral_style();
else if (!strcmp(command,"dimension")) dimension();
else if (!strcmp(command,"dump")) dump();
else if (!strcmp(command,"dump_modify")) dump_modify();
else if (!strcmp(command,"fix")) fix();
else if (!strcmp(command,"fix_modify")) fix_modify();
else if (!strcmp(command,"group")) group_command();
else if (!strcmp(command,"improper_coeff")) improper_coeff();
else if (!strcmp(command,"improper_style")) improper_style();
else if (!strcmp(command,"kspace_modify")) kspace_modify();
else if (!strcmp(command,"kspace_style")) kspace_style();
else if (!strcmp(command,"lattice")) lattice();
else if (!strcmp(command,"mass")) mass();
else if (!strcmp(command,"min_modify")) min_modify();
else if (!strcmp(command,"min_style")) min_style();
else if (!strcmp(command,"neigh_modify")) neigh_modify();
else if (!strcmp(command,"neighbor")) neighbor_command();
else if (!strcmp(command,"newton")) newton();
else if (!strcmp(command,"package")) package();
else if (!strcmp(command,"pair_coeff")) pair_coeff();
else if (!strcmp(command,"pair_modify")) pair_modify();
else if (!strcmp(command,"pair_style")) pair_style();
else if (!strcmp(command,"pair_write")) pair_write();
else if (!strcmp(command,"processors")) processors();
else if (!strcmp(command,"region")) region();
else if (!strcmp(command,"reset_timestep")) reset_timestep();
else if (!strcmp(command,"restart")) restart();
else if (!strcmp(command,"run_style")) run_style();
else if (!strcmp(command,"special_bonds")) special_bonds();
else if (!strcmp(command,"suffix")) suffix();
else if (!strcmp(command,"thermo")) thermo();
else if (!strcmp(command,"thermo_modify")) thermo_modify();
else if (!strcmp(command,"thermo_style")) thermo_style();
else if (!strcmp(command,"timestep")) timestep();
+ else if (!strcmp(command,"timers")) timers();
else if (!strcmp(command,"uncompute")) uncompute();
else if (!strcmp(command,"undump")) undump();
else if (!strcmp(command,"unfix")) unfix();
else if (!strcmp(command,"units")) units();
else flag = 0;
// return if command was listed above
if (flag) return 0;
// invoke commands added via style_command.h
if (command_map->find(command) != command_map->end()) {
CommandCreator command_creator = (*command_map)[command];
command_creator(lmp,narg,arg);
return 0;
}
// unrecognized command
return -1;
}
/* ----------------------------------------------------------------------
one instance per command in style_command.h
------------------------------------------------------------------------- */
template <typename T>
void Input::command_creator(LAMMPS *lmp, int narg, char **arg)
{
T cmd(lmp);
cmd.command(narg,arg);
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void Input::clear()
{
if (narg > 0) error->all(FLERR,"Illegal clear command");
lmp->destroy();
lmp->create();
lmp->post_create();
}
/* ---------------------------------------------------------------------- */
void Input::echo()
{
if (narg != 1) error->all(FLERR,"Illegal echo command");
if (strcmp(arg[0],"none") == 0) {
echo_screen = 0;
echo_log = 0;
} else if (strcmp(arg[0],"screen") == 0) {
echo_screen = 1;
echo_log = 0;
} else if (strcmp(arg[0],"log") == 0) {
echo_screen = 0;
echo_log = 1;
} else if (strcmp(arg[0],"both") == 0) {
echo_screen = 1;
echo_log = 1;
} else error->all(FLERR,"Illegal echo command");
}
/* ---------------------------------------------------------------------- */
void Input::ifthenelse()
{
if (narg < 3) error->all(FLERR,"Illegal if command");
// substitute for variables in Boolean expression for "if"
// in case expression was enclosed in quotes
// must substitute on copy of arg else will step on subsequent args
int n = strlen(arg[0]) + 1;
if (n > maxline) reallocate(line,maxline,n);
strcpy(line,arg[0]);
substitute(line,work,maxline,maxwork,0);
// evaluate Boolean expression for "if"
double btest = variable->evaluate_boolean(line);
// bound "then" commands
if (strcmp(arg[1],"then") != 0) error->all(FLERR,"Illegal if command");
int first = 2;
int iarg = first;
while (iarg < narg &&
(strcmp(arg[iarg],"elif") != 0 && strcmp(arg[iarg],"else") != 0))
iarg++;
int last = iarg-1;
// execute "then" commands
// make copies of all arg string commands
// required because re-parsing a command via one() will wipe out args
if (btest != 0.0) {
int ncommands = last-first + 1;
if (ncommands <= 0) error->all(FLERR,"Illegal if command");
char **commands = new char*[ncommands];
ncommands = 0;
for (int i = first; i <= last; i++) {
int n = strlen(arg[i]) + 1;
if (n == 1) error->all(FLERR,"Illegal if command");
commands[ncommands] = new char[n];
strcpy(commands[ncommands],arg[i]);
ncommands++;
}
ifthenelse_flag = 1;
for (int i = 0; i < ncommands; i++) one(commands[i]);
ifthenelse_flag = 0;
for (int i = 0; i < ncommands; i++) delete [] commands[i];
delete [] commands;
return;
}
// done if no "elif" or "else"
if (iarg == narg) return;
// check "elif" or "else" until find commands to execute
// substitute for variables and evaluate Boolean expression for "elif"
// must substitute on copy of arg else will step on subsequent args
// bound and execute "elif" or "else" commands
while (1) {
if (iarg+2 > narg) error->all(FLERR,"Illegal if command");
if (strcmp(arg[iarg],"elif") == 0) {
n = strlen(arg[iarg+1]) + 1;
if (n > maxline) reallocate(line,maxline,n);
strcpy(line,arg[iarg+1]);
substitute(line,work,maxline,maxwork,0);
btest = variable->evaluate_boolean(line);
first = iarg+2;
} else {
btest = 1.0;
first = iarg+1;
}
iarg = first;
while (iarg < narg &&
(strcmp(arg[iarg],"elif") != 0 && strcmp(arg[iarg],"else") != 0))
iarg++;
last = iarg-1;
if (btest == 0.0) continue;
int ncommands = last-first + 1;
if (ncommands <= 0) error->all(FLERR,"Illegal if command");
char **commands = new char*[ncommands];
ncommands = 0;
for (int i = first; i <= last; i++) {
int n = strlen(arg[i]) + 1;
if (n == 1) error->all(FLERR,"Illegal if command");
commands[ncommands] = new char[n];
strcpy(commands[ncommands],arg[i]);
ncommands++;
}
// execute the list of commands
ifthenelse_flag = 1;
for (int i = 0; i < ncommands; i++) one(commands[i]);
ifthenelse_flag = 0;
// clean up
for (int i = 0; i < ncommands; i++) delete [] commands[i];
delete [] commands;
return;
}
}
/* ---------------------------------------------------------------------- */
void Input::include()
{
if (narg != 1) error->all(FLERR,"Illegal include command");
// do not allow include inside an if command
// NOTE: this check will fail if a 2nd if command was inside the if command
// and came before the include
if (ifthenelse_flag)
error->all(FLERR,"Cannot use include command within an if command");
if (me == 0) {
if (nfile == maxfile) {
maxfile++;
infiles = (FILE **)
memory->srealloc(infiles,maxfile*sizeof(FILE *),"input:infiles");
}
infile = fopen(arg[0],"r");
if (infile == NULL) {
char str[128];
sprintf(str,"Cannot open input script %s",arg[0]);
error->one(FLERR,str);
}
infiles[nfile++] = infile;
}
}
/* ---------------------------------------------------------------------- */
void Input::jump()
{
if (narg < 1 || narg > 2) error->all(FLERR,"Illegal jump command");
if (jump_skip) {
jump_skip = 0;
return;
}
if (me == 0) {
if (strcmp(arg[0],"SELF") == 0) rewind(infile);
else {
if (infile != stdin) fclose(infile);
infile = fopen(arg[0],"r");
if (infile == NULL) {
char str[128];
sprintf(str,"Cannot open input script %s",arg[0]);
error->one(FLERR,str);
}
infiles[nfile-1] = infile;
}
}
if (narg == 2) {
label_active = 1;
if (labelstr) delete [] labelstr;
int n = strlen(arg[1]) + 1;
labelstr = new char[n];
strcpy(labelstr,arg[1]);
}
}
/* ---------------------------------------------------------------------- */
void Input::label()
{
if (narg != 1) error->all(FLERR,"Illegal label command");
if (label_active && strcmp(labelstr,arg[0]) == 0) label_active = 0;
}
/* ---------------------------------------------------------------------- */
void Input::log()
{
if (narg > 2) error->all(FLERR,"Illegal log command");
int appendflag = 0;
if (narg == 2) {
if (strcmp(arg[1],"append") == 0) appendflag = 1;
else error->all(FLERR,"Illegal log command");
}
if (me == 0) {
if (logfile) fclose(logfile);
if (strcmp(arg[0],"none") == 0) logfile = NULL;
else {
if (appendflag) logfile = fopen(arg[0],"a");
else logfile = fopen(arg[0],"w");
if (logfile == NULL) {
char str[128];
sprintf(str,"Cannot open logfile %s",arg[0]);
error->one(FLERR,str);
}
}
if (universe->nworlds == 1) universe->ulogfile = logfile;
}
}
/* ---------------------------------------------------------------------- */
void Input::next_command()
{
if (variable->next(narg,arg)) jump_skip = 1;
}
/* ---------------------------------------------------------------------- */
void Input::partition()
{
if (narg < 3) error->all(FLERR,"Illegal partition command");
int yesflag;
if (strcmp(arg[0],"yes") == 0) yesflag = 1;
else if (strcmp(arg[0],"no") == 0) yesflag = 0;
else error->all(FLERR,"Illegal partition command");
int ilo,ihi;
force->bounds(arg[1],universe->nworlds,ilo,ihi);
// copy original line to copy, since will use strtok() on it
// ptr = start of 4th word
strcpy(copy,line);
char *ptr = strtok(copy," \t\n\r\f");
ptr = strtok(NULL," \t\n\r\f");
ptr = strtok(NULL," \t\n\r\f");
ptr += strlen(ptr) + 1;
ptr += strspn(ptr," \t\n\r\f");
// execute the remaining command line on requested partitions
if (yesflag) {
if (universe->iworld+1 >= ilo && universe->iworld+1 <= ihi) one(ptr);
} else {
if (universe->iworld+1 < ilo || universe->iworld+1 > ihi) one(ptr);
}
}
/* ---------------------------------------------------------------------- */
void Input::print()
{
if (narg < 1) error->all(FLERR,"Illegal print command");
// copy 1st arg back into line (copy is being used)
// check maxline since arg[0] could have been exanded by variables
// substitute for $ variables (no printing) and print arg
int n = strlen(arg[0]) + 1;
if (n > maxline) reallocate(line,maxline,n);
strcpy(line,arg[0]);
substitute(line,work,maxline,maxwork,0);
// parse optional args
FILE *fp = NULL;
int screenflag = 1;
int iarg = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"file") == 0 || strcmp(arg[iarg],"append") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal print command");
if (me == 0) {
if (strcmp(arg[iarg],"file") == 0) fp = fopen(arg[iarg+1],"w");
else fp = fopen(arg[iarg+1],"a");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open print file %s",arg[iarg+1]);
error->one(FLERR,str);
}
}
iarg += 2;
} else if (strcmp(arg[iarg],"screen") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal print command");
if (strcmp(arg[iarg+1],"yes") == 0) screenflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) screenflag = 0;
else error->all(FLERR,"Illegal print command");
iarg += 2;
} else error->all(FLERR,"Illegal print command");
}
if (me == 0) {
if (screenflag && screen) fprintf(screen,"%s\n",line);
if (screenflag && logfile) fprintf(logfile,"%s\n",line);
if (fp) {
fprintf(fp,"%s\n",line);
fclose(fp);
}
}
}
/* ---------------------------------------------------------------------- */
void Input::quit()
{
if (narg) error->all(FLERR,"Illegal quit command");
error->done();
}
/* ---------------------------------------------------------------------- */
void Input::shell()
{
if (narg < 1) error->all(FLERR,"Illegal shell command");
if (strcmp(arg[0],"cd") == 0) {
- if (narg != 2) error->all(FLERR,"Illegal shell command");
+ if (narg != 2) error->all(FLERR,"Illegal shell cd command");
chdir(arg[1]);
} else if (strcmp(arg[0],"mkdir") == 0) {
- if (narg < 2) error->all(FLERR,"Illegal shell command");
+ if (narg < 2) error->all(FLERR,"Illegal shell mkdir command");
#if !defined(WINDOWS) && !defined(__MINGW32__)
if (me == 0)
for (int i = 1; i < narg; i++)
mkdir(arg[i], S_IRWXU | S_IRGRP | S_IXGRP);
#endif
} else if (strcmp(arg[0],"mv") == 0) {
- if (narg != 3) error->all(FLERR,"Illegal shell command");
+ if (narg != 3) error->all(FLERR,"Illegal shell mv command");
if (me == 0) rename(arg[1],arg[2]);
} else if (strcmp(arg[0],"rm") == 0) {
- if (narg < 2) error->all(FLERR,"Illegal shell command");
+ if (narg < 2) error->all(FLERR,"Illegal shell rm command");
if (me == 0)
for (int i = 1; i < narg; i++) unlink(arg[i]);
} else if (strcmp(arg[0],"rmdir") == 0) {
- if (narg < 2) error->all(FLERR,"Illegal shell command");
+ if (narg < 2) error->all(FLERR,"Illegal shell rmdir command");
if (me == 0)
for (int i = 1; i < narg; i++) rmdir(arg[i]);
+ } else if (strcmp(arg[0],"setenv") == 0) {
+ if (narg != 3) error->all(FLERR,"Illegal shell setenv command");
+ setenv(arg[1],arg[2],1);
+
// use work string to concat args back into one string separated by spaces
// invoke string in shell via system()
} else {
int n = 0;
for (int i = 0; i < narg; i++) n += strlen(arg[i]) + 1;
if (n > maxwork) reallocate(work,maxwork,n);
strcpy(work,arg[0]);
for (int i = 1; i < narg; i++) {
strcat(work," ");
strcat(work,arg[i]);
}
if (me == 0) system(work);
}
}
/* ---------------------------------------------------------------------- */
void Input::variable_command()
{
variable->set(narg,arg);
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
one function for each LAMMPS-specific input script command
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void Input::angle_coeff()
{
if (domain->box_exist == 0)
error->all(FLERR,"Angle_coeff command before simulation box is defined");
if (force->angle == NULL)
error->all(FLERR,"Angle_coeff command before angle_style is defined");
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Angle_coeff command when no angles allowed");
force->angle->coeff(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::angle_style()
{
if (narg < 1) error->all(FLERR,"Illegal angle_style command");
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Angle_style command when no angles allowed");
force->create_angle(arg[0],lmp->suffix);
if (force->angle) force->angle->settings(narg-1,&arg[1]);
}
/* ---------------------------------------------------------------------- */
void Input::atom_modify()
{
atom->modify_params(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::atom_style()
{
if (narg < 1) error->all(FLERR,"Illegal atom_style command");
if (domain->box_exist)
error->all(FLERR,"Atom_style command after simulation box is defined");
atom->create_avec(arg[0],narg-1,&arg[1],lmp->suffix);
}
/* ---------------------------------------------------------------------- */
void Input::bond_coeff()
{
if (domain->box_exist == 0)
error->all(FLERR,"Bond_coeff command before simulation box is defined");
if (force->bond == NULL)
error->all(FLERR,"Bond_coeff command before bond_style is defined");
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Bond_coeff command when no bonds allowed");
force->bond->coeff(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::bond_style()
{
if (narg < 1) error->all(FLERR,"Illegal bond_style command");
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Bond_style command when no bonds allowed");
force->create_bond(arg[0],lmp->suffix);
if (force->bond) force->bond->settings(narg-1,&arg[1]);
}
/* ---------------------------------------------------------------------- */
void Input::boundary()
{
if (domain->box_exist)
error->all(FLERR,"Boundary command after simulation box is defined");
domain->set_boundary(narg,arg,0);
}
/* ---------------------------------------------------------------------- */
void Input::box()
{
if (domain->box_exist)
error->all(FLERR,"Box command after simulation box is defined");
domain->set_box(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::communicate()
{
comm->set(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::compute()
{
modify->add_compute(narg,arg,lmp->suffix);
}
/* ---------------------------------------------------------------------- */
void Input::compute_modify()
{
modify->modify_compute(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::dielectric()
{
if (narg != 1) error->all(FLERR,"Illegal dielectric command");
force->dielectric = force->numeric(FLERR,arg[0]);
}
/* ---------------------------------------------------------------------- */
void Input::dihedral_coeff()
{
if (domain->box_exist == 0)
error->all(FLERR,"Dihedral_coeff command before simulation box is defined");
if (force->dihedral == NULL)
error->all(FLERR,"Dihedral_coeff command before dihedral_style is defined");
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Dihedral_coeff command when no dihedrals allowed");
force->dihedral->coeff(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::dihedral_style()
{
if (narg < 1) error->all(FLERR,"Illegal dihedral_style command");
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Dihedral_style command when no dihedrals allowed");
force->create_dihedral(arg[0],lmp->suffix);
if (force->dihedral) force->dihedral->settings(narg-1,&arg[1]);
}
/* ---------------------------------------------------------------------- */
void Input::dimension()
{
if (narg != 1) error->all(FLERR,"Illegal dimension command");
if (domain->box_exist)
error->all(FLERR,"Dimension command after simulation box is defined");
domain->dimension = force->inumeric(FLERR,arg[0]);
if (domain->dimension != 2 && domain->dimension != 3)
error->all(FLERR,"Illegal dimension command");
// must reset default extra_dof of all computes
// since some were created before dimension command is encountered
for (int i = 0; i < modify->ncompute; i++)
modify->compute[i]->reset_extra_dof();
}
/* ---------------------------------------------------------------------- */
void Input::dump()
{
output->add_dump(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::dump_modify()
{
output->modify_dump(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::fix()
{
modify->add_fix(narg,arg,lmp->suffix);
}
/* ---------------------------------------------------------------------- */
void Input::fix_modify()
{
modify->modify_fix(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::group_command()
{
group->assign(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::improper_coeff()
{
if (domain->box_exist == 0)
error->all(FLERR,"Improper_coeff command before simulation box is defined");
if (force->improper == NULL)
error->all(FLERR,"Improper_coeff command before improper_style is defined");
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Improper_coeff command when no impropers allowed");
force->improper->coeff(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::improper_style()
{
if (narg < 1) error->all(FLERR,"Illegal improper_style command");
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Improper_style command when no impropers allowed");
force->create_improper(arg[0],lmp->suffix);
if (force->improper) force->improper->settings(narg-1,&arg[1]);
}
/* ---------------------------------------------------------------------- */
void Input::kspace_modify()
{
if (force->kspace == NULL)
error->all(FLERR,"KSpace style has not yet been set");
force->kspace->modify_params(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::kspace_style()
{
force->create_kspace(narg,arg,lmp->suffix);
}
/* ---------------------------------------------------------------------- */
void Input::lattice()
{
domain->set_lattice(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::mass()
{
if (narg != 2) error->all(FLERR,"Illegal mass command");
if (domain->box_exist == 0)
error->all(FLERR,"Mass command before simulation box is defined");
atom->set_mass(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::min_modify()
{
update->minimize->modify_params(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::min_style()
{
if (domain->box_exist == 0)
error->all(FLERR,"Min_style command before simulation box is defined");
update->create_minimize(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::neigh_modify()
{
neighbor->modify_params(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::neighbor_command()
{
neighbor->set(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::newton()
{
int newton_pair,newton_bond;
if (narg == 1) {
if (strcmp(arg[0],"off") == 0) newton_pair = newton_bond = 0;
else if (strcmp(arg[0],"on") == 0) newton_pair = newton_bond = 1;
else error->all(FLERR,"Illegal newton command");
} else if (narg == 2) {
if (strcmp(arg[0],"off") == 0) newton_pair = 0;
else if (strcmp(arg[0],"on") == 0) newton_pair= 1;
else error->all(FLERR,"Illegal newton command");
if (strcmp(arg[1],"off") == 0) newton_bond = 0;
else if (strcmp(arg[1],"on") == 0) newton_bond = 1;
else error->all(FLERR,"Illegal newton command");
} else error->all(FLERR,"Illegal newton command");
force->newton_pair = newton_pair;
if (newton_bond == 0) {
if (domain->box_exist && force->newton_bond == 1)
error->all(FLERR,"Newton bond change after simulation box is defined");
force->newton_bond = 0;
} else {
if (domain->box_exist && force->newton_bond == 0)
error->all(FLERR,"Newton bond change after simulation box is defined");
force->newton_bond = 1;
}
if (newton_pair || newton_bond) force->newton = 1;
else force->newton = 0;
}
/* ---------------------------------------------------------------------- */
void Input::package()
{
if (domain->box_exist)
error->all(FLERR,"Package command after simulation box is defined");
if (narg < 1) error->all(FLERR,"Illegal package command");
if (strcmp(arg[0],"cuda") == 0) {
if (!lmp->cuda)
error->all(FLERR,"Package cuda command without USER-CUDA installed");
lmp->cuda->accelerator(narg-1,&arg[1]);
} else if (strcmp(arg[0],"gpu") == 0) {
char **fixarg = new char*[2+narg];
fixarg[0] = (char *) "package_gpu";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "GPU";
for (int i = 1; i < narg; i++) fixarg[i+2] = arg[i];
modify->add_fix(2+narg,fixarg,NULL);
delete [] fixarg;
force->newton_pair = 0;
} else if (strcmp(arg[0],"omp") == 0) {
char **fixarg = new char*[2+narg];
fixarg[0] = (char *) "package_omp";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "OMP";
for (int i = 1; i < narg; i++) fixarg[i+2] = arg[i];
modify->add_fix(2+narg,fixarg,NULL);
delete [] fixarg;
} else error->all(FLERR,"Illegal package command");
}
/* ---------------------------------------------------------------------- */
void Input::pair_coeff()
{
if (domain->box_exist == 0)
error->all(FLERR,"Pair_coeff command before simulation box is defined");
if (force->pair == NULL)
error->all(FLERR,"Pair_coeff command before pair_style is defined");
force->pair->coeff(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::pair_modify()
{
if (force->pair == NULL)
error->all(FLERR,"Pair_modify command before pair_style is defined");
force->pair->modify_params(narg,arg);
}
/* ----------------------------------------------------------------------
if old pair style exists and new style is same, just change settings
else create new pair class
------------------------------------------------------------------------- */
void Input::pair_style()
{
if (narg < 1) error->all(FLERR,"Illegal pair_style command");
if (force->pair && strcmp(arg[0],force->pair_style) == 0) {
force->pair->settings(narg-1,&arg[1]);
return;
}
force->create_pair(arg[0],lmp->suffix);
if (force->pair) force->pair->settings(narg-1,&arg[1]);
}
/* ---------------------------------------------------------------------- */
void Input::pair_write()
{
if (force->pair == NULL)
error->all(FLERR,"Pair_write command before pair_style is defined");
force->pair->write_file(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::processors()
{
if (domain->box_exist)
error->all(FLERR,"Processors command after simulation box is defined");
comm->set_processors(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::region()
{
domain->add_region(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::reset_timestep()
{
update->reset_timestep(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::restart()
{
output->create_restart(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::run_style()
{
if (domain->box_exist == 0)
error->all(FLERR,"Run_style command before simulation box is defined");
update->create_integrate(narg,arg,lmp->suffix);
}
/* ---------------------------------------------------------------------- */
void Input::special_bonds()
{
// store 1-3,1-4 and dihedral/extra flag values before change
// change in 1-2 coeffs will not change the special list
double lj2 = force->special_lj[2];
double lj3 = force->special_lj[3];
double coul2 = force->special_coul[2];
double coul3 = force->special_coul[3];
int angle = force->special_angle;
int dihedral = force->special_dihedral;
int extra = force->special_extra;
force->set_special(narg,arg);
// if simulation box defined and saved values changed, redo special list
if (domain->box_exist && atom->molecular) {
if (lj2 != force->special_lj[2] || lj3 != force->special_lj[3] ||
coul2 != force->special_coul[2] || coul3 != force->special_coul[3] ||
angle != force->special_angle ||
dihedral != force->special_dihedral ||
extra != force->special_extra) {
Special special(lmp);
special.build();
}
}
}
/* ---------------------------------------------------------------------- */
void Input::suffix()
{
if (narg != 1) error->all(FLERR,"Illegal suffix command");
if (strcmp(arg[0],"off") == 0) lmp->suffix_enable = 0;
else if (strcmp(arg[0],"on") == 0) lmp->suffix_enable = 1;
else {
delete [] lmp->suffix;
int n = strlen(arg[0]) + 1;
lmp->suffix = new char[n];
strcpy(lmp->suffix,arg[0]);
lmp->suffix_enable = 1;
}
}
/* ---------------------------------------------------------------------- */
void Input::thermo()
{
output->set_thermo(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::thermo_modify()
{
output->thermo->modify_params(narg,arg);
}
/* ---------------------------------------------------------------------- */
void Input::thermo_style()
{
output->create_thermo(narg,arg);
}
/* ---------------------------------------------------------------------- */
+void Input::timers()
+{
+ timer->modify_params(narg,arg);
+}
+
+/* ---------------------------------------------------------------------- */
+
void Input::timestep()
{
if (narg != 1) error->all(FLERR,"Illegal timestep command");
update->dt = force->numeric(FLERR,arg[0]);
}
/* ---------------------------------------------------------------------- */
void Input::uncompute()
{
if (narg != 1) error->all(FLERR,"Illegal uncompute command");
modify->delete_compute(arg[0]);
}
/* ---------------------------------------------------------------------- */
void Input::undump()
{
if (narg != 1) error->all(FLERR,"Illegal undump command");
output->delete_dump(arg[0]);
}
/* ---------------------------------------------------------------------- */
void Input::unfix()
{
if (narg != 1) error->all(FLERR,"Illegal unfix command");
modify->delete_fix(arg[0]);
}
/* ---------------------------------------------------------------------- */
void Input::units()
{
if (narg != 1) error->all(FLERR,"Illegal units command");
if (domain->box_exist)
error->all(FLERR,"Units command after simulation box is defined");
update->set_units(arg[0]);
}