jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #603 from jrgissing/Ncreate_atoms 

add exactly N particles to available lattice points
This commit is contained in:
Axel Kohlmeyer 2020-01-22 15:36:54 -05:00 committed by GitHub
parent 64a4bf6a42 ac143dbdb7
commit 9d06430894
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
16 changed files with 496 additions and 406 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Errors_messages.rst
View File

@ -286,6 +286,9 @@ Doc page with :doc:`WARNING messages <Errors_warnings>`
*Attempting to rescale a 0.0 temperature* *Attempting to rescale a 0.0 temperature*
Cannot rescale a temperature that is already 0.0. Cannot rescale a temperature that is already 0.0.
*Attempting to insert more particles than available lattice points*
Self-explanatory.
*Bad FENE bond* *Bad FENE bond*
Two atoms in a FENE bond have become so far apart that the bond cannot Two atoms in a FENE bond have become so far apart that the bond cannot
be computed. be computed.

3
doc/src/Errors_warnings.rst
View File

@ -683,6 +683,9 @@ This will most likely cause errors in kinetic fluctuations.
*Slab correction not needed for MSM* *Slab correction not needed for MSM*
Slab correction is intended to be used with Ewald or PPPM and is not needed by MSM. Slab correction is intended to be used with Ewald or PPPM and is not needed by MSM.
*Specifying an 'subset' value of '0' is equivalent to no 'subset' keyword*
Self-explanatory.
*System is not charge neutral, net charge = %g* *System is not charge neutral, net charge = %g*
The total charge on all atoms on the system is not 0.0. The total charge on all atoms on the system is not 0.0.
For some KSpace solvers this is only a warning. For some KSpace solvers this is only a warning.

22
doc/src/create_atoms.rst
View File

@ -27,7 +27,7 @@ Syntax
region-ID = create atoms within this region, use NULL for entire simulation box region-ID = create atoms within this region, use NULL for entire simulation box
* zero or more keyword/value pairs may be appended * zero or more keyword/value pairs may be appended
* keyword = *mol* or *basis* or *remap* or *var* or *set* or *units* * keyword = *mol* or *basis* or *ratio* or *subset* or *remap* or *var* or *set* or *rotate* or *units*
.. parsed-literal:: .. parsed-literal::
@ -37,6 +37,12 @@ Syntax
*basis* values = M itype *basis* values = M itype
M = which basis atom M = which basis atom
itype = atom type (1-N) to assign to this basis atom itype = atom type (1-N) to assign to this basis atom
*ratio* values = frac seed
frac = fraction of lattice sites (0 to 1) to populate randomly
seed = random # seed (positive integer)
*subset* values = Nsubset seed
Nsubset = # of lattice sites to populate randomly
seed = random # seed (positive integer)
*remap* value = *yes* or *no* *remap* value = *yes* or *no*
*var* value = name = variable name to evaluate for test of atom creation *var* value = name = variable name to evaluate for test of atom creation
*set* values = dim name *set* values = dim name
@ -59,6 +65,7 @@ Examples
create_atoms 1 box create_atoms 1 box
create_atoms 3 region regsphere basis 2 3 create_atoms 3 region regsphere basis 2 3
create_atoms 3 region regsphere basis 2 3 ratio 0.5 74637
create_atoms 3 single 0 0 5 create_atoms 3 single 0 0 5
create_atoms 1 box var v set x xpos set y ypos create_atoms 1 box var v set x xpos set y ypos
@ -214,6 +221,19 @@ command for specifics on how basis atoms are defined for the unit cell
of the lattice. By default, all created atoms are assigned the of the lattice. By default, all created atoms are assigned the
argument *type* as their atom type. argument *type* as their atom type.
The *ratio* and *subset* keywords can be used in conjunction with the
*box* or *region* styles to limit the total number of particles
inserted. The lattice defines a set of *Nlatt* eligible sites for
inserting particles, which may be limited by the *region* style or the
*var* and *set* keywords. For the *ratio* keyword only the specified
fraction of them (0 <= *frac* <= 1) will be assigned particles. For
the *subset* keyword only the specified *Nsubset* of them will be
assigned particles. In both cases the assigned lattice sites are
chosen randomly. An iterative algorithm is used which insures the
correct number of particles are inserted, in a perfectly random
fashion. Which lattice sites are selected will change with the number
of processors used.
The *remap* keyword only applies to the *single* style. If it is set The *remap* keyword only applies to the *single* style. If it is set
to *yes*\ , then if the specified position is outside the simulation to *yes*\ , then if the specified position is outside the simulation
box, it will mapped back into the box, assuming the relevant box, it will mapped back into the box, assuming the relevant

68
doc/src/set.rst
View File

@ -14,7 +14,7 @@ Syntax
* style = *atom* or *type* or *mol* or *group* or *region* * style = *atom* or *type* or *mol* or *group* or *region*
* ID = atom ID range or type range or mol ID range or group ID or region ID * ID = atom ID range or type range or mol ID range or group ID or region ID
* one or more keyword/value pairs may be appended * one or more keyword/value pairs may be appended
* keyword = *type* or *type/fraction* or *mol* or *x* or *y* or *z* or *charge* or *dipole* or *dipole/random* or *quat* or *spin* or *spin/random* or *quat* or *quat/random* or *diameter* or *shape* or *length* or *tri* or *theta* or *theta/random* or *angmom* or *omega* or *mass* or *density* or *density/disc* or *volume* or *image* or *bond* or *angle* or *dihedral* or *improper* or *meso/e* or *meso/cv* or *meso/rho* or *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or *edpd/temp* or *edpd/cv* or *cc* or *i\_name* or *d\_name* * keyword = *type* or *type/fraction* or *type/ratio* or *type/subset* or *mol* or *x* or *y* or *z* or *charge* or *dipole* or *dipole/random* or *quat* or *spin* or *spin/random* or *quat* or *quat/random* or *diameter* or *shape* or *length* or *tri* or *theta* or *theta/random* or *angmom* or *omega* or *mass* or *density* or *density/disc* or *volume* or *image* or *bond* or *angle* or *dihedral* or *improper* or *meso/e* or *meso/cv* or *meso/rho* or *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or *edpd/temp* or *edpd/cv* or *cc* or *i\_name* or *d\_name*
.. parsed-literal:: .. parsed-literal::
@ -22,7 +22,15 @@ Syntax
value can be an atom-style variable (see below) value can be an atom-style variable (see below)
*type/fraction* values = type fraction seed *type/fraction* values = type fraction seed
type = new atom type type = new atom type
fraction = fraction of selected atoms to set to new atom type fraction = approximate fraction of selected atoms to set to new atom type
seed = random # seed (positive integer)
*type/ratio* values = type fraction seed
type = new atom type
fraction = exact fraction of selected atoms to set to new atom type
seed = random # seed (positive integer)
*type/subset* values = type Nsubset seed
type = new atom type
Nsubset = exact number of selected atoms to set to new atom type
seed = random # seed (positive integer) seed = random # seed (positive integer)
*mol* value = molecule ID *mol* value = molecule ID
value can be an atom-style variable (see below) value can be an atom-style variable (see below)
@ -184,15 +192,16 @@ This section describes the keyword options for which properties to
change, for the selected atoms. change, for the selected atoms.
Note that except where explicitly prohibited below, all of the Note that except where explicitly prohibited below, all of the
keywords allow an :doc:`atom-style or atomfile-style variable <variable>` to be used as the specified value(s). If the keywords allow an :doc:`atom-style or atomfile-style variable
value is a variable, it should be specified as v\_name, where name is <variable>` to be used as the specified value(s). If the value is a
the variable name. In this case, the variable will be evaluated, and variable, it should be specified as v\_name, where name is the
its resulting per-atom value used to determine the value assigned to variable name. In this case, the variable will be evaluated, and its
each selected atom. Note that the per-atom value from the variable resulting per-atom value used to determine the value assigned to each
will be ignored for atoms that are not selected via the *style* and selected atom. Note that the per-atom value from the variable will be
*ID* settings explained above. A simple way to use per-atom values ignored for atoms that are not selected via the *style* and *ID*
from the variable to reset a property for all atoms is to use style settings explained above. A simple way to use per-atom values from
*atom* with *ID* = "\*"; this selects all atom IDs. the variable to reset a property for all atoms is to use style *atom*
with *ID* = "\*"; this selects all atom IDs.
Atom-style variables can specify formulas with various mathematical Atom-style variables can specify formulas with various mathematical
functions, and include :doc:`thermo\_style <thermo_style>` command functions, and include :doc:`thermo\_style <thermo_style>` command
@ -220,23 +229,36 @@ command.
Keyword *type/fraction* sets the atom type for a fraction of the Keyword *type/fraction* sets the atom type for a fraction of the
selected atoms. The actual number of atoms changed is not guaranteed selected atoms. The actual number of atoms changed is not guaranteed
to be exactly the requested fraction, but should be statistically to be exactly the specified fraction (0 <= *fraction* <= 1), but
close. Random numbers are used in such a way that a particular atom should be statistically close. Random numbers are used in such a way
is changed or not changed, regardless of how many processors are being that a particular atom is changed or not changed, regardless of how
used. This keyword does not allow use of an atom-style variable. many processors are being used. This keyword does not allow use of an
atom-style variable.
Keyword *mol* sets the molecule ID for all selected atoms. The :doc:`atom style <atom_style>` being used must support the use of molecule Keywords *type/ratio* and *type/subset" also set the atom type for a
IDs. fraction of the selected atoms. The actual number of atoms changed
will be exactly the requested number. For *type/ratio* the specified
fraction (0 <= *fraction* <= 1) determines the number. For
*type/subset*, the specified *Nsubset* is the number. An iterative
algorithm is used which insures the correct number of atoms are
selected, in a perfectly random fashion. Which atoms are selected
will change with the number of processors used. These keywords do not
allow use of an atom-style variable.
Keywords *x*\ , *y*\ , *z*\ , and *charge* set the coordinates or charge of Keyword *mol* sets the molecule ID for all selected atoms. The
all selected atoms. For *charge*\ , the :doc:`atom style <atom_style>` :doc:`atom style <atom_style>` being used must support the use of
being used must support the use of atomic charge. Keywords *vx*\ , *vy*\ , molecule IDs.
and *vz* set the velocities of all selected atoms.
Keywords *x*\ , *y*\ , *z*\ , and *charge* set the coordinates or
charge of all selected atoms. For *charge*\ , the :doc:`atom style
<atom_style>` being used must support the use of atomic
charge. Keywords *vx*\ , *vy*\ , and *vz* set the velocities of all
selected atoms.
Keyword *dipole* uses the specified x,y,z values as components of a Keyword *dipole* uses the specified x,y,z values as components of a
vector to set as the orientation of the dipole moment vectors of the vector to set as the orientation of the dipole moment vectors of the
selected atoms. The magnitude of the dipole moment is set selected atoms. The magnitude of the dipole moment is set by the
by the length of this orientation vector. length of this orientation vector.
Keyword *dipole/random* randomizes the orientation of the dipole Keyword *dipole/random* randomizes the orientation of the dipole
moment vectors for the selected atoms and sets the magnitude of each moment vectors for the selected atoms and sets the magnitude of each

334
doc/txt/create_atoms.txt
View File

@ -1,334 +0,0 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Commands_all.html)
:line
create_atoms command :h3
[Syntax:]
create_atoms type style args keyword values ... :pre
type = atom type (1-Ntypes) of atoms to create (offset for molecule creation) :ulb,l
style = {box} or {region} or {single} or {random} :l
{box} args = none
{region} args = region-ID
region-ID = particles will only be created if contained in the region
{single} args = x y z
x,y,z = coordinates of a single particle (distance units)
{random} args = N seed region-ID
N = number of particles to create
seed = random # seed (positive integer)
region-ID = create atoms within this region, use NULL for entire simulation box :pre
zero or more keyword/value pairs may be appended :l
keyword = {mol} or {basis} or {remap} or {var} or {set} or {units} :l
{mol} value = template-ID seed
template-ID = ID of molecule template specified in a separate "molecule"_molecule.html command
seed = random # seed (positive integer)
{basis} values = M itype
M = which basis atom
itype = atom type (1-N) to assign to this basis atom
{remap} value = {yes} or {no}
{var} value = name = variable name to evaluate for test of atom creation
{set} values = dim name
dim = {x} or {y} or {z}
name = name of variable to set with x, y, or z atom position
{rotate} values = theta Rx Ry Rz
theta = rotation angle for single molecule (degrees)
Rx,Ry,Rz = rotation vector for single molecule
{units} value = {lattice} or {box}
{lattice} = the geometry is defined in lattice units
{box} = the geometry is defined in simulation box units :pre
:ule
[Examples:]
create_atoms 1 box
create_atoms 3 region regsphere basis 2 3
create_atoms 3 single 0 0 5
create_atoms 1 box var v set x xpos set y ypos :pre
[Description:]
This command creates atoms (or molecules) on a lattice, or a single
atom (or molecule), or a random collection of atoms (or molecules), as
an alternative to reading in their coordinates explicitly via a
"read_data"_read_data.html or "read_restart"_read_restart.html
command. A simulation box must already exist, which is typically
created via the "create_box"_create_box.html command. Before using
this command, a lattice must also be defined using the
"lattice"_lattice.html command, unless you specify the {single} style
with units = box or the {random} style. For the remainder of this doc
page, a created atom or molecule is referred to as a "particle".
If created particles are individual atoms, they are assigned the
specified atom {type}, though this can be altered via the {basis}
keyword as discussed below. If molecules are being created, the type
of each atom in the created molecule is specified in the file read by
the "molecule"_molecule.html command, and those values are added to
the specified atom {type}. E.g. if {type} = 2, and the file specifies
atom types 1,2,3, then each created molecule will have atom types
3,4,5.
For the {box} style, the create_atoms command fills the entire
simulation box with particles on the lattice. If your simulation box
is periodic, you should insure its size is a multiple of the lattice
spacings, to avoid unwanted atom overlaps at the box boundaries. If
your box is periodic and a multiple of the lattice spacing in a
particular dimension, LAMMPS is careful to put exactly one particle at
the boundary (on either side of the box), not zero or two.
For the {region} style, a geometric volume is filled with particles on
the lattice. This volume what is inside the simulation box and is
also consistent with the region volume. See the "region"_region.html
command for details. Note that a region can be specified so that its
"volume" is either inside or outside a geometric boundary. Also note
that if your region is the same size as a periodic simulation box (in
some dimension), LAMMPS does not implement the same logic described
above as for the {box} style, to insure exactly one particle at
periodic boundaries. if this is what you desire, you should either
use the {box} style, or tweak the region size to get precisely the
particles you want.
For the {single} style, a single particle is added to the system at
the specified coordinates. This can be useful for debugging purposes
or to create a tiny system with a handful of particles at specified
positions.
For the {random} style, N particles are added to the system at
randomly generated coordinates, which can be useful for generating an
amorphous system. The particles are created one by one using the
specified random number {seed}, resulting in the same set of particles
coordinates, independent of how many processors are being used in the
simulation. If the {region-ID} argument is specified as NULL, then
the created particles will be anywhere in the simulation box. If a
{region-ID} is specified, a geometric volume is filled which is both
inside the simulation box and is also consistent with the region
volume. See the "region"_region.html command for details. Note that
a region can be specified so that its "volume" is either inside or
outside a geometric boundary.
NOTE: Particles generated by the {random} style will typically be
highly overlapped which will cause many interatomic potentials to
compute large energies and forces. Thus you should either perform an
"energy minimization"_minimize.html or run dynamics with "fix
nve/limit"_fix_nve_limit.html to equilibrate such a system, before
running normal dynamics.
Note that this command adds particles to those that already exist.
This means it can be used to add particles to a system previously read
in from a data or restart file. Or the create_atoms command can be
used multiple times, to add multiple sets of particles to the
simulation. For example, grain boundaries can be created, by
interleaving create_atoms with "lattice"_lattice.html commands
specifying different orientations. By using the create_atoms command
in conjunction with the "delete_atoms"_delete_atoms.html command,
reasonably complex geometries can be created, or a protein can be
solvated with a surrounding box of water molecules.
In all these cases, care should be taken to insure that new atoms do
not overlap existing atoms inappropriately, especially if molecules
are being added. The "delete_atoms"_delete_atoms.html command can be
used to remove overlapping atoms or molecules.
NOTE: You cannot use any of the styles explained above to create atoms
that are outside the simulation box; they will just be ignored by
LAMMPS. This is true even if you are using shrink-wrapped box
boundaries, as specified by the "boundary"_boundary.html command.
However, you can first use the "change_box"_change_box.html command to
temporarily expand the box, then add atoms via create_atoms, then
finally use change_box command again if needed to re-shrink-wrap the
new atoms. See the "change_box"_change_box.html doc page for an
example of how to do this, using the create_atoms {single} style to
insert a new atom outside the current simulation box.
:line
Individual atoms are inserted by this command, unless the {mol}
keyword is used. It specifies a {template-ID} previously defined
using the "molecule"_molecule.html command, which reads a file that
defines the molecule. The coordinates, atom types, charges, etc, as
well as any bond/angle/etc and special neighbor information for the
molecule can be specified in the molecule file. See the
"molecule"_molecule.html command for details. The only settings
required to be in this file are the coordinates and types of atoms in
the molecule.
Using a lattice to add molecules, e.g. via the {box} or {region} or
{single} styles, is exactly the same as adding atoms on lattice
points, except that entire molecules are added at each point, i.e. on
the point defined by each basis atom in the unit cell as it tiles the
simulation box or region. This is done by placing the geometric
center of the molecule at the lattice point, and giving the molecule a
random orientation about the point. The random {seed} specified with
the {mol} keyword is used for this operation, and the random numbers
generated by each processor are different. This means the coordinates
of individual atoms (in the molecules) will be different when running
on different numbers of processors, unlike when atoms are being
created in parallel.
Also note that because of the random rotations, it may be important to
use a lattice with a large enough spacing that adjacent molecules will
not overlap, regardless of their relative orientations.
NOTE: If the "create_box"_create_box.html command is used to create
the simulation box, followed by the create_atoms command with its
{mol} option for adding molecules, then you typically need to use the
optional keywords allowed by the "create_box"_create_box.html command
for extra bonds (angles,etc) or extra special neighbors. This is
because by default, the "create_box"_create_box.html command sets up a
non-molecular system which doesn't allow molecules to be added.
:line
This is the meaning of the other allowed keywords.
The {basis} keyword is only used when atoms (not molecules) are being
created. It specifies an atom type that will be assigned to specific
basis atoms as they are created. See the "lattice"_lattice.html
command for specifics on how basis atoms are defined for the unit cell
of the lattice. By default, all created atoms are assigned the
argument {type} as their atom type.
The {remap} keyword only applies to the {single} style. If it is set
to {yes}, then if the specified position is outside the simulation
box, it will mapped back into the box, assuming the relevant
dimensions are periodic. If it is set to {no}, no remapping is done
and no particle is created if its position is outside the box.
The {var} and {set} keywords can be used together to provide a
criterion for accepting or rejecting the addition of an individual
atom, based on its coordinates. The {name} specified for the {var}
keyword is the name of an "equal-style variable"_variable.html which
should evaluate to a zero or non-zero value based on one or two or
three variables which will store the x, y, or z coordinates of an atom
(one variable per coordinate). If used, these other variables must be
"internal-style variables"_variable.html defined in the input script;
their initial numeric value can be anything. They must be
internal-style variables, because this command resets their values
directly. The {set} keyword is used to identify the names of these
other variables, one variable for the x-coordinate of a created atom,
one for y, and one for z.
When an atom is created, its x,y,z coordinates become the values for
any {set} variable that is defined. The {var} variable is then
evaluated. If the returned value is 0.0, the atom is not created. If
it is non-zero, the atom is created.
As an example, these commands can be used in a 2d simulation, to
create a sinusoidal surface. Note that the surface is "rough" due to
individual lattice points being "above" or "below" the mathematical
expression for the sinusoidal curve. If a finer lattice were used,
the sinusoid would appear to be "smoother". Also note the use of the
"xlat" and "ylat" "thermo_style"_thermo_style.html keywords which
converts lattice spacings to distance. Click on the image for a
larger version.
dimension 2
variable x equal 100
variable y equal 25
lattice hex 0.8442
region box block 0 $x 0 $y -0.5 0.5
create_box 1 box :pre
variable xx internal 0.0
variable yy internal 0.0
variable v equal "(0.2*v_y*ylat * cos(v_xx/xlat * 2.0*PI*4.0/v_x) + 0.5*v_y*ylat - v_yy) > 0.0"
create_atoms 1 box var v set x xx set y yy
write_dump all atom sinusoid.lammpstrj :pre
:c,image(JPG/sinusoid_small.jpg,JPG/sinusoid.jpg)
The {rotate} keyword allows specification of the orientation
at which molecules are inserted. The axis of rotation is
determined by the rotation vector (Rx,Ry,Rz) that goes through the
insertion point. The specified {theta} determines the angle of
rotation around that axis. Note that the direction of rotation for
the atoms around the rotation axis is consistent with the right-hand
rule: if your right-hand's thumb points along {R}, then your fingers
wrap around the axis in the direction of rotation.
The {units} keyword determines the meaning of the distance units used
to specify the coordinates of the one particle created by the {single}
style. A {box} value selects standard distance units as defined by
the "units"_units.html command, e.g. Angstroms for units = real or
metal. A {lattice} value means the distance units are in lattice
spacings.
:line
Atom IDs are assigned to created atoms in the following way. The
collection of created atoms are assigned consecutive IDs that start
immediately following the largest atom ID existing before the
create_atoms command was invoked. This is done by the processor's
communicating the number of atoms they each own, the first processor
numbering its atoms from 1 to N1, the second processor from N1+1 to
N2, etc. Where N1 = number of atoms owned by the first processor, N2
= number owned by the second processor, etc. Thus when the same
simulation is performed on different numbers of processors, there is
no guarantee a particular created atom will be assigned the same ID in
both simulations. If molecules are being created, molecule IDs are
assigned to created molecules in a similar fashion.
Aside from their ID, atom type, and xyz position, other properties of
created atoms are set to default values, depending on which quantities
are defined by the chosen "atom style"_atom_style.html. See the "atom
style"_atom_style.html command for more details. See the
"set"_set.html and "velocity"_velocity.html commands for info on how
to change these values.
charge = 0.0
dipole moment magnitude = 0.0
diameter = 1.0
shape = 0.0 0.0 0.0
density = 1.0
volume = 1.0
velocity = 0.0 0.0 0.0
angular velocity = 0.0 0.0 0.0
angular momentum = 0.0 0.0 0.0
quaternion = (1,0,0,0)
bonds, angles, dihedrals, impropers = none :ul
If molecules are being created, these defaults can be overridden by
values specified in the file read by the "molecule"_molecule.html
command. E.g. the file typically defines bonds (angles,etc) between
atoms in the molecule, and can optionally define charges on each atom.
Note that the {sphere} atom style sets the default particle diameter
to 1.0 as well as the density. This means the mass for the particle
is not 1.0, but is PI/6 * diameter^3 = 0.5236.
Note that the {ellipsoid} atom style sets the default particle shape
to (0.0 0.0 0.0) and the density to 1.0 which means it is a point
particle, not an ellipsoid, and has a mass of 1.0.
Note that the {peri} style sets the default volume and density to 1.0
and thus also set the mass for the particle to 1.0.
The "set"_set.html command can be used to override many of these
default settings.
:line
[Restrictions:]
An "atom_style"_atom_style.html must be previously defined to use this
command.
A rotation vector specified for a single molecule must be in
the z-direction for a 2d model.
[Related commands:]
"lattice"_lattice.html, "region"_region.html, "create_box"_create_box.html,
"read_data"_read_data.html, "read_restart"_read_restart.html
[Default:]
The default for the {basis} keyword is that all created atoms are
assigned the argument {type} as their atom type (when single atoms are
being created). The other defaults are {remap} = no, {rotate} =
random, and {units} = lattice.

15
src/atom_vec.cpp
View File

@ -19,6 +19,8 @@
#include "error.h" #include "error.h"
#include "utils.h" #include "utils.h"
#include "comm.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
#define DELTA 16384 #define DELTA 16384
@ -108,6 +110,19 @@ int AtomVec::grow_nmax_bonus(int nmax_bonus)
return nmax_bonus; return nmax_bonus;
} }
/* ----------------------------------------------------------------------
roundup N so it is a multiple of DELTA
error if N exceeds 32-bit int, since will be used as arg to grow()
------------------------------------------------------------------------- */
bigint AtomVec::roundup(bigint n)
{
if (n % DELTA) n = n/DELTA * DELTA + DELTA;
if (n > MAXSMALLINT)
error->one(FLERR,"Too many atoms created on one or more procs");
return n;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
unpack one line from Velocities section of data file unpack one line from Velocities section of data file
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */

1
src/atom_vec.h
View File

@ -57,6 +57,7 @@ class AtomVec : protected Pointers {
virtual void grow(int) = 0; virtual void grow(int) = 0;
virtual void grow_reset() = 0; virtual void grow_reset() = 0;
bigint roundup(bigint);
virtual void copy(int, int, int) = 0; virtual void copy(int, int, int) = 0;
virtual void clear_bonus() {} virtual void clear_bonus() {}
virtual void force_clear(int, size_t) {} virtual void force_clear(int, size_t) {}

134
src/create_atoms.cpp
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory. See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author (ratio and subset) : Jake Gissinger (U Colorado)
------------------------------------------------------------------------- */
#include "create_atoms.h" #include "create_atoms.h"
#include <mpi.h> #include <mpi.h>
#include <cstring> #include <cstring>
@ -32,15 +36,19 @@
#include "math_extra.h" #include "math_extra.h"
#include "math_const.h" #include "math_const.h"
#include "error.h" #include "error.h"
#include "memory.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
#define BIG 1.0e30 #define BIG 1.0e30
#define EPSILON 1.0e-6 #define EPSILON 1.0e-6
#define LB_FACTOR 1.1
enum{BOX,REGION,SINGLE,RANDOM}; enum{BOX,REGION,SINGLE,RANDOM};
enum{ATOM,MOLECULE}; enum{ATOM,MOLECULE};
enum{COUNT,INSERT,INSERT_SELECTED};
enum{NONE,RATIO,SUBSET};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -50,6 +58,9 @@ CreateAtoms::CreateAtoms(LAMMPS *lmp) : Pointers(lmp) {}
void CreateAtoms::command(int narg, char **arg) void CreateAtoms::command(int narg, char **arg)
{ {
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
if (domain->box_exist == 0) if (domain->box_exist == 0)
error->all(FLERR,"Create_atoms command before simulation box is defined"); error->all(FLERR,"Create_atoms command before simulation box is defined");
if (modify->nfix_restart_peratom) if (modify->nfix_restart_peratom)
@ -107,6 +118,8 @@ void CreateAtoms::command(int narg, char **arg)
varflag = 0; varflag = 0;
vstr = xstr = ystr = zstr = NULL; vstr = xstr = ystr = zstr = NULL;
quatone[0] = quatone[1] = quatone[2] = 0.0; quatone[0] = quatone[1] = quatone[2] = 0.0;
subsetflag = NONE;
int subsetseed;
nbasis = domain->lattice->nbasis; nbasis = domain->lattice->nbasis;
basistype = new int[nbasis]; basistype = new int[nbasis];
@ -132,7 +145,7 @@ void CreateAtoms::command(int narg, char **arg)
int imol = atom->find_molecule(arg[iarg+1]); int imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1) error->all(FLERR,"Molecule template ID for " if (imol == -1) error->all(FLERR,"Molecule template ID for "
"create_atoms does not exist"); "create_atoms does not exist");
if (atom->molecules[imol]->nset > 1 && comm->me == 0) if (atom->molecules[imol]->nset > 1 && me == 0)
error->warning(FLERR,"Molecule template for " error->warning(FLERR,"Molecule template for "
"create_atoms has multiple molecules"); "create_atoms has multiple molecules");
mode = MOLECULE; mode = MOLECULE;
@ -187,6 +200,22 @@ void CreateAtoms::command(int narg, char **arg)
MathExtra::norm3(axisone); MathExtra::norm3(axisone);
MathExtra::axisangle_to_quat(axisone,thetaone,quatone); MathExtra::axisangle_to_quat(axisone,thetaone,quatone);
iarg += 5; iarg += 5;
} else if (strcmp(arg[iarg],"ratio") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command");
subsetflag = RATIO;
subsetfrac = force->numeric(FLERR,arg[iarg+1]);
subsetseed = force->inumeric(FLERR,arg[iarg+2]);
if (subsetfrac <= 0.0 || subsetfrac > 1.0 || subsetseed <= 0)
error->all(FLERR,"Illegal create_atoms command");
iarg += 3;
} else if (strcmp(arg[iarg],"subset") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command");
subsetflag = SUBSET;
nsubset = force->bnumeric(FLERR,arg[iarg+1]);
subsetseed = force->inumeric(FLERR,arg[iarg+2]);
if (nsubset <= 0 || subsetseed <= 0)
error->all(FLERR,"Illegal create_atoms command");
iarg += 3;
} else error->all(FLERR,"Illegal create_atoms command"); } else error->all(FLERR,"Illegal create_atoms command");
} }
@ -221,9 +250,12 @@ void CreateAtoms::command(int narg, char **arg)
// molecule random number generator, different for each proc // molecule random number generator, different for each proc
ranmol = new RanMars(lmp,molseed+comm->me); ranmol = new RanMars(lmp,molseed+me);
} }
ranlatt = NULL;
if (subsetflag != NONE) ranlatt = new RanMars(lmp,subsetseed+me);
// error check and further setup for variable test // error check and further setup for variable test
if (!vstr && (xstr || ystr || zstr)) if (!vstr && (xstr || ystr || zstr))
@ -512,6 +544,8 @@ void CreateAtoms::command(int narg, char **arg)
// clean up // clean up
delete ranmol; delete ranmol;
delete ranlatt;
if (domain->lattice) delete [] basistype; if (domain->lattice) delete [] basistype;
delete [] vstr; delete [] vstr;
delete [] xstr; delete [] xstr;
@ -536,7 +570,7 @@ void CreateAtoms::command(int narg, char **arg)
MPI_Barrier(world); MPI_Barrier(world);
double time2 = MPI_Wtime(); double time2 = MPI_Wtime();
if (comm->me == 0) { if (me == 0) {
if (screen) { if (screen) {
fprintf(screen,"Created " BIGINT_FORMAT " atoms\n", fprintf(screen,"Created " BIGINT_FORMAT " atoms\n",
atom->natoms-natoms_previous); atom->natoms-natoms_previous);
@ -757,7 +791,6 @@ void CreateAtoms::add_lattice()
// which can lead to missing atoms in rare cases // which can lead to missing atoms in rare cases
// extra decrement of lo if min < 0, since static_cast(-1.5) = -1 // extra decrement of lo if min < 0, since static_cast(-1.5) = -1
int ilo,ihi,jlo,jhi,klo,khi;
ilo = static_cast<int> (xmin) - 1; ilo = static_cast<int> (xmin) - 1;
jlo = static_cast<int> (ymin) - 1; jlo = static_cast<int> (ymin) - 1;
klo = static_cast<int> (zmin) - 1; klo = static_cast<int> (zmin) - 1;
@ -769,14 +802,73 @@ void CreateAtoms::add_lattice()
if (ymin < 0.0) jlo--; if (ymin < 0.0) jlo--;
if (zmin < 0.0) klo--; if (zmin < 0.0) klo--;
// iterate on 3d periodic lattice of unit cells using loop bounds // count lattice sites on each proc
// iterate on nbasis atoms in each unit cell
// convert lattice coords to box coords nlatt_overflow = 0;
// add atom or molecule (on each basis point) if it meets all criteria loop_lattice(COUNT);
// nadd = # of atoms each proc will insert (estimated if subsetflag)
int overflow;
MPI_Allreduce(&nlatt_overflow,&overflow,1,MPI_INT,MPI_SUM,world);
if (overflow)
error->all(FLERR,"Create_atoms lattice size overflow on 1 or more procs");
bigint nadd;
if (subsetflag == NONE) {
if (nprocs == 1) nadd = nlatt;
else nadd = static_cast<bigint> (LB_FACTOR * nlatt);
} else {
bigint bnlatt = nlatt;
bigint bnlattall;
MPI_Allreduce(&bnlatt,&bnlattall,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (subsetflag == RATIO)
nsubset = static_cast<bigint> (subsetfrac * bnlattall);
if (nsubset > bnlattall)
error->all(FLERR,"Create_atoms subset size > # of lattice sites");
if (nprocs == 1) nadd = nsubset;
else nadd = static_cast<bigint> (LB_FACTOR * nsubset/bnlattall * nlatt);
}
// allocate atom arrays to size N, rounded up by AtomVec->DELTA
bigint nbig = atom->avec->roundup(nadd + atom->nlocal);
int n = static_cast<int> (nbig);
atom->avec->grow(n);
// add atoms or molecules
// if no subset: add to all lattice sites
// if subset: count lattice sites, select random subset, then add
if (subsetflag == NONE) loop_lattice(INSERT);
else {
memory->create(flag,nlatt,"create_atoms:flag");
memory->create(next,nlatt,"create_atoms:next");
ranlatt->select_subset(nsubset,nlatt,flag,next);
loop_lattice(INSERT_SELECTED);
memory->destroy(flag);
memory->destroy(next);
}
}
/* ----------------------------------------------------------------------
iterate on 3d periodic lattice of unit cells using loop bounds
iterate on nbasis atoms in each unit cell
convert lattice coords to box coords
check if lattice point meets all criteria to be added
perform action on atom or molecule (on each basis point) if meets all criteria
actions = add, count, add if flagged
------------------------------------------------------------------------- */
void CreateAtoms::loop_lattice(int action)
{
int i,j,k,m;
const double * const * const basis = domain->lattice->basis; const double * const * const basis = domain->lattice->basis;
int i,j,k,m; nlatt = 0;
for (k = klo; k <= khi; k++) { for (k = klo; k <= khi; k++) {
for (j = jlo; j <= jhi; j++) { for (j = jlo; j <= jhi; j++) {
for (i = ilo; i <= ihi; i++) { for (i = ilo; i <= ihi; i++) {
@ -812,19 +904,31 @@ void CreateAtoms::add_lattice()
coord[1] < sublo[1] || coord[1] >= subhi[1] || coord[1] < sublo[1] || coord[1] >= subhi[1] ||
coord[2] < sublo[2] || coord[2] >= subhi[2]) continue; coord[2] < sublo[2] || coord[2] >= subhi[2]) continue;
// add the atom or entire molecule to my list of atoms // this proc owns the lattice site
// perform action: add, just count, add if flagged
// add = add an atom or entire molecule to my list of atoms
if (mode == ATOM) atom->avec->create_atom(basistype[m],x); if (action == INSERT) {
else if (quatone[0] == 0 && quatone[1] == 0 && quatone[2] == 0) if (mode == ATOM) atom->avec->create_atom(basistype[m],x);
add_molecule(x); else if (quatone[0] == 0 && quatone[1] == 0 && quatone[2] == 0)
else add_molecule(x,quatone); add_molecule(x);
else add_molecule(x,quatone);
} else if (action == COUNT) {
if (nlatt == MAXSMALLINT) nlatt_overflow = 1;
} else if (action == INSERT_SELECTED && flag[nlatt]) {
if (mode == ATOM) atom->avec->create_atom(basistype[m],x);
else if (quatone[0] == 0 && quatone[1] == 0 && quatone[2] == 0)
add_molecule(x);
else add_molecule(x,quatone);
}
nlatt++;
} }
} }
} }
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
add a randomly rotated molecule with its center at center add a randomly rotated molecule with its center at center
if quat_user set, perform requested rotation if quat_user set, perform requested rotation

26
src/create_atoms.h
View File

@ -30,17 +30,30 @@ class CreateAtoms : protected Pointers {
void command(int, char **); void command(int, char **);
private: private:
int me,nprocs;
int ntype,style,mode,nregion,nbasis,nrandom,seed; int ntype,style,mode,nregion,nbasis,nrandom,seed;
int remapflag;
int subsetflag;
bigint nsubset;
double subsetfrac;
int *basistype; int *basistype;
double xone[3],quatone[4]; double xone[3],quatone[4];
int remapflag;
int varflag,vvar,xvar,yvar,zvar; int varflag,vvar,xvar,yvar,zvar;
char *vstr,*xstr,*ystr,*zstr; char *vstr,*xstr,*ystr,*zstr;
char *xstr_copy,*ystr_copy,*zstr_copy; char *xstr_copy,*ystr_copy,*zstr_copy;
int ilo,ihi,jlo,jhi,klo,khi;
int nlatt; // number of owned lattice sites
int nlatt_overflow; // 1 if local nlatt exceeds a 32-bit int
int *flag; // flag subset of particles to insert on lattice
int *next;
class Molecule *onemol; class Molecule *onemol;
class RanMars *ranmol; class RanMars *ranmol;
class RanMars *ranlatt;
int triclinic; int triclinic;
double sublo[3],subhi[3]; // epsilon-extended proc sub-box for adding atoms double sublo[3],subhi[3]; // epsilon-extended proc sub-box for adding atoms
@ -48,8 +61,9 @@ class CreateAtoms : protected Pointers {
void add_single(); void add_single();
void add_random(); void add_random();
void add_lattice(); void add_lattice();
void loop_lattice(int);
void add_molecule(double *, double * = NULL); void add_molecule(double *, double * = NULL);
int vartest(double *); // evaluate a variable with new atom position int vartest(double *); // evaluate a variable with new atom position
}; };
} }
@ -152,4 +166,12 @@ E: No overlap of box and region for create_atoms
Self-explanatory. Self-explanatory.
E: Attempting to insert more particles than available lattice points
Self-explanatory.
W: Specifying an 'subset' value of '0' is equivalent to no 'subset' keyword
Self-explanatory.
*/ */

132
src/random_mars.cpp
View File

@ -21,6 +21,8 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{ADD,SUBTRACT};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp), RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp),
@ -149,7 +151,7 @@ double RanMars::besselexp(double theta, double alpha, double cp)
{ {
double first,v1,v2; double first,v1,v2;
if (theta < 0.0 || alpha < 0.0 || alpha > 1.0) if (theta < 0.0 || alpha < 0.0 || alpha < 1.0)
error->all(FLERR,"Invalid Bessel exponential distribution parameters"); error->all(FLERR,"Invalid Bessel exponential distribution parameters");
v1 = uniform(); v1 = uniform();
@ -166,3 +168,131 @@ double RanMars::besselexp(double theta, double alpha, double cp)
return first; return first;
} }
/* ----------------------------------------------------------------------
select random subset of size Ntarget out of Ntotal items
Ntotal = sum of Nmine across all procs
mark,next = vectors of length Nmine
return mark = 0 for unselected item, 1 for selected item
next = work vector used to store linked lists for 2 active sets of items
IMPORTANT: this method must be called simultaneously by all procs
------------------------------------------------------------------------- */
void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next)
{
int mode,index,oldindex,newvalue,nflip,which,niter;
int active[2],first[2],last[2];
int newactive[2],newfirst[2],newlast[2];
bigint nmark,nactive,nactiveall,nflipall,bnflip;
bigint activeall[2],bsum[3],bsumall[3];
double thresh;
active[0] = nmine;
active[1] = 0;
first[0] = 0;
first[1] = -1;
last[0] = nmine-1;
last[1] = -1;
bigint bnmine = nmine;
bigint bnall;
MPI_Allreduce(&bnmine,&bnall,1,MPI_LMP_BIGINT,MPI_SUM,world);
activeall[0] = bnall;
for (int i = 0; i < nmine; i++) mark[i] = 0;
for (int i = 0; i < nmine; i++) next[i] = i+1;
next[nmine-1] = -1;
nmark = 0;
niter = 0;
while (nmark != ntarget) {
// choose to ADD or SUBTRACT from current nmark
// thresh = desired flips / size of active set
// nactive = size of current active set, only for debug output below
if (ntarget-nmark > 0) {
mode = ADD;
// nactive = active[mode];
thresh = 1.0 * (ntarget-nmark) / activeall[mode];
} else {
mode = SUBTRACT;
// nactive = active[mode];
thresh = 1.0 * (nmark-ntarget) / activeall[mode];
}
// bound thresh for RNG accuracy
thresh = MAX(thresh,0.01);
thresh = MIN(thresh,0.99);
// new empty active sets for next iteration
newactive[0] = newactive[1] = 0;
newfirst[0] = newfirst[1] = -1;
newlast[0] = newlast[1] = -1;
// index = first value in ADD or SUBTRACT set
if (mode == ADD) newvalue = 1;
else if (mode == SUBTRACT) newvalue = 0;
index = first[mode];
// flip marks from 0 -> 1 (ADD) or 1 -> 0 (SUBTRACT)
// loop over active set via next vector = linked list
// flip each value based on RN < thresh
nflip = 0;
while (index >= 0) {
if (uniform() < thresh) {
mark[index] = newvalue;
nflip++;
}
oldindex = index;
index = next[index];
// oldindex can now be appended to a new active set
// which = which of two new active sets to append to
which = mark[oldindex];
newactive[which]++;
if (newfirst[which] < 0) newfirst[which] = oldindex;
if (newlast[which] >= 0) next[newlast[which]] = oldindex;
newlast[which] = oldindex;
next[oldindex] = -1;
// set active sets for next iteration to the new ones
// next vector is already updated
active[0] = newactive[0];
active[1] = newactive[1];
first[0] = newfirst[0];
first[1] = newfirst[1];
last[0] = newlast[0];
last[1] = newlast[1];
}
// update nmark and activeall
bsum[0] = nflip;
bsum[1] = active[0];
bsum[2] = active[1];
bsum[3] = nactive;
MPI_Allreduce(&bsum,&bsumall,4,MPI_LMP_BIGINT,MPI_SUM,world);
nflipall = bsumall[0];
activeall[0] = bsumall[1];
activeall[1] = bsumall[2];
nactiveall = bsumall[3];
if (mode == ADD) nmark += nflipall;
else if (mode == SUBTRACT) nmark -= nflipall;
niter++;
// DEBUG output of stats
//if (comm->me == 0) printf("%d %ld %ld %g %ld\n",
// niter,nmark,nactiveall,thresh,nflipall);
}
}

1
src/random_mars.h
View File

@ -27,6 +27,7 @@ class RanMars : protected Pointers {
double gaussian(double mu, double sigma); double gaussian(double mu, double sigma);
double rayleigh(double sigma); double rayleigh(double sigma);
double besselexp(double theta, double alpha, double cp); double besselexp(double theta, double alpha, double cp);
void select_subset(bigint, int, int *, int *);
private: private:
int save; int save;

6
src/read_data.cpp
View File

@ -440,6 +440,12 @@ void ReadData::command(int narg, char **arg)
atom->allocate_type_arrays(); atom->allocate_type_arrays();
atom->deallocate_topology(); atom->deallocate_topology();
// allocate atom arrays to N, rounded up by AtomVec->DELTA
bigint nbig = n;
nbig = atom->avec->roundup(nbig);
n = static_cast<int> (nbig);
atom->avec->grow(n); atom->avec->grow(n);
domain->boxlo[0] = boxlo[0]; domain->boxhi[0] = boxhi[0]; domain->boxlo[0] = boxlo[0]; domain->boxhi[0] = boxhi[0];

12
src/read_restart.cpp
View File

@ -165,6 +165,12 @@ void ReadRestart::command(int narg, char **arg)
atom->allocate_type_arrays(); atom->allocate_type_arrays();
atom->deallocate_topology(); atom->deallocate_topology();
// allocate atom arrays to size N, rounded up by AtomVec->DELTA
bigint nbig = n;
nbig = atom->avec->roundup(nbig);
n = static_cast<int> (nbig);
atom->avec->grow(n); atom->avec->grow(n);
n = atom->nmax; n = atom->nmax;
@ -211,13 +217,17 @@ void ReadRestart::command(int narg, char **arg)
memory->create(buf,assignedChunkSize,"read_restart:buf"); memory->create(buf,assignedChunkSize,"read_restart:buf");
mpiio->read((headerOffset+assignedChunkOffset),assignedChunkSize,buf); mpiio->read((headerOffset+assignedChunkOffset),assignedChunkSize,buf);
mpiio->close(); mpiio->close();
if (!nextra) { // We can actually calculate number of atoms from assignedChunkSize
// can calculate number of atoms from assignedChunkSize
if (!nextra) {
atom->nlocal = 1; // temporarily claim there is one atom... atom->nlocal = 1; // temporarily claim there is one atom...
int perAtomSize = avec->size_restart(); // ...so we can get its size int perAtomSize = avec->size_restart(); // ...so we can get its size
atom->nlocal = 0; // restore nlocal to zero atoms atom->nlocal = 0; // restore nlocal to zero atoms
int atomCt = (int) (assignedChunkSize / perAtomSize); int atomCt = (int) (assignedChunkSize / perAtomSize);
if (atomCt > atom->nmax) avec->grow(atomCt); if (atomCt > atom->nmax) avec->grow(atomCt);
} }
m = 0; m = 0;
while (m < assignedChunkSize) m += avec->unpack_restart(&buf[m]); while (m < assignedChunkSize) m += avec->unpack_restart(&buf[m]);
} }

6
src/replicate.cpp
View File

@ -223,6 +223,12 @@ void Replicate::command(int narg, char **arg)
else n = static_cast<int> (LB_FACTOR * atom->natoms / nprocs); else n = static_cast<int> (LB_FACTOR * atom->natoms / nprocs);
atom->allocate_type_arrays(); atom->allocate_type_arrays();
// allocate atom arrays to size N, rounded up by AtomVec->DELTA
bigint nbig = n;
nbig = atom->avec->roundup(nbig);
n = static_cast<int> (nbig);
atom->avec->grow(n); atom->avec->grow(n);
n = atom->nmax; n = atom->nmax;

136
src/set.cpp
View File

@ -30,6 +30,7 @@
#include "input.h" #include "input.h"
#include "variable.h" #include "variable.h"
#include "random_park.h" #include "random_park.h"
#include "random_mars.h"
#include "math_extra.h" #include "math_extra.h"
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
@ -41,7 +42,8 @@ using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI, enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET,
MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM, DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM,
THETA,THETA_RANDOM,ANGMOM,OMEGA, THETA,THETA_RANDOM,ANGMOM,OMEGA,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
@ -109,6 +111,34 @@ void Set::command(int narg, char **arg)
setrandom(TYPE_FRACTION); setrandom(TYPE_FRACTION);
iarg += 4; iarg += 4;
} else if (strcmp(arg[iarg],"type/ratio") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
newtype = force->inumeric(FLERR,arg[iarg+1]);
fraction = force->numeric(FLERR,arg[iarg+2]);
ivalue = force->inumeric(FLERR,arg[iarg+3]);
if (newtype <= 0 || newtype > atom->ntypes)
error->all(FLERR,"Invalid value in set command");
if (fraction < 0.0 || fraction > 1.0)
error->all(FLERR,"Invalid value in set command");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
setrandom(TYPE_RATIO);
iarg += 4;
} else if (strcmp(arg[iarg],"type/subset") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
newtype = force->inumeric(FLERR,arg[iarg+1]);
nsubset = force->bnumeric(FLERR,arg[iarg+2]);
ivalue = force->inumeric(FLERR,arg[iarg+3]);
if (newtype <= 0 || newtype > atom->ntypes)
error->all(FLERR,"Invalid value in set command");
if (nsubset < 0)
error->all(FLERR,"Invalid value in set command");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
setrandom(TYPE_SUBSET);
iarg += 4;
} else if (strcmp(arg[iarg],"mol") == 0) { } else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
@ -1001,23 +1031,72 @@ void Set::setrandom(int keyword)
AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri"); AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body"); AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body");
RanPark *random = new RanPark(lmp,1);
double **x = atom->x; double **x = atom->x;
int seed = ivalue; int seed = ivalue;
// set fraction of atom types to newtype RanPark *ranpark = new RanPark(lmp,1);
RanMars *ranmars = new RanMars(lmp,seed + comm->me);
// set approx fraction of atom types to newtype
if (keyword == TYPE_FRACTION) { if (keyword == TYPE_FRACTION) {
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (select[i]) { if (select[i]) {
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
if (random->uniform() > fraction) continue; if (ranpark->uniform() > fraction) continue;
atom->type[i] = newtype; atom->type[i] = newtype;
count++; count++;
} }
// set exact count of atom types to newtype
// for TYPE_RATIO, exact = fraction out of total eligible
// for TYPE_SUBSET, exact = nsubset out of total eligible
} else if (keyword == TYPE_RATIO || keyword == TYPE_SUBSET) {
int nlocal = atom->nlocal;
// count = number of eligible atoms I own
count = 0;
for (i = 0; i < nlocal; i++)
if (select[i]) count++;
// convert specified fraction to nsubset
bigint bcount = count;
bigint allcount;
MPI_Allreduce(&bcount,&allcount,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (keyword == TYPE_RATIO) {
nsubset = static_cast<bigint> (fraction * allcount);
} else if (keyword == TYPE_SUBSET) {
if (nsubset > allcount)
error->all(FLERR,"Set type/subset value exceeds eligible atoms");
}
// make selection
int *flag = memory->create(flag,count,"set:flag");
int *work = memory->create(work,count,"set:work");
ranmars->select_subset(nsubset,count,flag,work);
// change types of selected atoms
count = 0;
for (i = 0; i < nlocal; i++)
if (select[i] && flag[i]) {
atom->type[i] = newtype;
count++;
}
// clean up
memory->destroy(flag);
memory->destroy(work);
// set dipole moments to random orientations in 3d or 2d // set dipole moments to random orientations in 3d or 2d
// dipole length is determined by dipole type array // dipole length is determined by dipole type array
@ -1030,10 +1109,10 @@ void Set::setrandom(int keyword)
if (domain->dimension == 3) { if (domain->dimension == 3) {
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (select[i]) { if (select[i]) {
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
mu[i][0] = random->uniform() - 0.5; mu[i][0] = ranpark->uniform() - 0.5;
mu[i][1] = random->uniform() - 0.5; mu[i][1] = ranpark->uniform() - 0.5;
mu[i][2] = random->uniform() - 0.5; mu[i][2] = ranpark->uniform() - 0.5;
msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
scale = dvalue/sqrt(msq); scale = dvalue/sqrt(msq);
mu[i][0] *= scale; mu[i][0] *= scale;
@ -1046,9 +1125,9 @@ void Set::setrandom(int keyword)
} else { } else {
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (select[i]) { if (select[i]) {
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
mu[i][0] = random->uniform() - 0.5; mu[i][0] = ranpark->uniform() - 0.5;
mu[i][1] = random->uniform() - 0.5; mu[i][1] = ranpark->uniform() - 0.5;
mu[i][2] = 0.0; mu[i][2] = 0.0;
msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1]; msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1];
scale = dvalue/sqrt(msq); scale = dvalue/sqrt(msq);
@ -1072,10 +1151,10 @@ void Set::setrandom(int keyword)
if (domain->dimension == 3) { if (domain->dimension == 3) {
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (select[i]) { if (select[i]) {
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
sp[i][0] = random->uniform() - 0.5; sp[i][0] = ranpark->uniform() - 0.5;
sp[i][1] = random->uniform() - 0.5; sp[i][1] = ranpark->uniform() - 0.5;
sp[i][2] = random->uniform() - 0.5; sp[i][2] = ranpark->uniform() - 0.5;
sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1] + sp[i][2]*sp[i][2]; sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1] + sp[i][2]*sp[i][2];
scale = 1.0/sqrt(sp_sq); scale = 1.0/sqrt(sp_sq);
sp[i][0] *= scale; sp[i][0] *= scale;
@ -1088,9 +1167,9 @@ void Set::setrandom(int keyword)
} else { } else {
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (select[i]) { if (select[i]) {
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
sp[i][0] = random->uniform() - 0.5; sp[i][0] = ranpark->uniform() - 0.5;
sp[i][1] = random->uniform() - 0.5; sp[i][1] = ranpark->uniform() - 0.5;
sp[i][2] = 0.0; sp[i][2] = 0.0;
sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1]; sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1];
scale = 1.0/sqrt(sp_sq); scale = 1.0/sqrt(sp_sq);
@ -1120,12 +1199,12 @@ void Set::setrandom(int keyword)
else else
error->one(FLERR,"Cannot set quaternion for atom that has none"); error->one(FLERR,"Cannot set quaternion for atom that has none");
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
s = random->uniform(); s = ranpark->uniform();
t1 = sqrt(1.0-s); t1 = sqrt(1.0-s);
t2 = sqrt(s); t2 = sqrt(s);
theta1 = 2.0*MY_PI*random->uniform(); theta1 = 2.0*MY_PI*ranpark->uniform();
theta2 = 2.0*MY_PI*random->uniform(); theta2 = 2.0*MY_PI*ranpark->uniform();
quat[0] = cos(theta2)*t2; quat[0] = cos(theta2)*t2;
quat[1] = sin(theta1)*t1; quat[1] = sin(theta1)*t1;
quat[2] = cos(theta1)*t1; quat[2] = cos(theta1)*t1;
@ -1144,8 +1223,8 @@ void Set::setrandom(int keyword)
else else
error->one(FLERR,"Cannot set quaternion for atom that has none"); error->one(FLERR,"Cannot set quaternion for atom that has none");
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
theta2 = MY_PI*random->uniform(); theta2 = MY_PI*ranpark->uniform();
quat[0] = cos(theta2); quat[0] = cos(theta2);
quat[1] = 0.0; quat[1] = 0.0;
quat[2] = 0.0; quat[2] = 0.0;
@ -1162,14 +1241,15 @@ void Set::setrandom(int keyword)
if (select[i]) { if (select[i]) {
if (atom->line[i] < 0) if (atom->line[i] < 0)
error->one(FLERR,"Cannot set theta for atom that is not a line"); error->one(FLERR,"Cannot set theta for atom that is not a line");
random->reset(seed,x[i]); ranpark->reset(seed,x[i]);
avec_line->bonus[atom->line[i]].theta = MY_2PI*random->uniform(); avec_line->bonus[atom->line[i]].theta = MY_2PI*ranpark->uniform();
count++; count++;
} }
} }
} }
delete random; delete ranpark;
delete ranmars;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

3
src/set.h
View File

@ -34,8 +34,9 @@ class Set : protected Pointers {
int *select; int *select;
int style,ivalue,newtype,count,index_custom; int style,ivalue,newtype,count,index_custom;
int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag; int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag;
double dvalue,xvalue,yvalue,zvalue,wvalue,fraction;
int cc_index; int cc_index;
bigint nsubset;
double dvalue,xvalue,yvalue,zvalue,wvalue,fraction;
int varflag,varflag1,varflag2,varflag3,varflag4; int varflag,varflag1,varflag2,varflag3,varflag4;
int ivar1,ivar2,ivar3,ivar4; int ivar1,ivar2,ivar3,ivar4;