Merge pull request #857 from lammps/resetids

new reset_ids command and dump_modify delay option
This commit is contained in:
Steve Plimpton
2018-03-30 09:45:02 -06:00
committed by GitHub
16 changed files with 438 additions and 42 deletions

View File

@ -497,6 +497,7 @@ in the command's documentation.
"region"_region.html,
"replicate"_replicate.html,
"rerun"_rerun.html,
"reset_ids"_reset_ids.html,
"reset_timestep"_reset_timestep.html,
"restart"_restart.html,
"run"_run.html,

View File

@ -82,6 +82,7 @@ Commands :h1
region
replicate
rerun
reset_ids
reset_timestep
restart
run

View File

@ -15,8 +15,9 @@ compute ID group-ID displace/atom :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
displace/atom = style name of this compute command :l
zero or more keyword/arg pairs may be appended :l
keyword = {refresh} :l
{refresh} arg = per-atom variable ID :pre
keyword = {refresh} :
{replace} arg = name of per-atom variable :pre
:ule
[Examples:]
@ -75,36 +76,41 @@ a distance {Dhop}. For any snapshot we only want to output atoms that
have hopped since the last snapshot. This can be accomplished with
something like the following commands:
write_dump all custom tmp.dump id type x y z # see comment below :pre
variable Dhop equal 0.6
variable check atom "c_dsp[4] > v_Dhop"
compute dsp all displace/atom refresh check
dump 1 all custom 20 tmp.dump id type x y z
dump_modify 1 append yes thresh c_dsp[4] > ${Dhop} refresh c_dsp :pre
dump 1 all custom 100 tmp.dump id type x y z
dump_modify 1 append yes thresh c_dsp[4] > ${Dhop} &
refresh c_dsp delay 100 :pre
The "dump_modify thresh"_dump_modify.html command will cause only
atoms that have displaced more than 0.6 Angstroms to be output on a
given snapshot (assuming metal units). The dump_modify {refresh}
option triggers a call to this compute at the end of every dump.
The "dump_modify thresh"_dump_modify.html command will only ouptut
atoms that have displaced more than 0.6 Angstroms on each snapshot
(assuming metal units). The dump_modify {refresh} option triggers a
call to this compute at the end of every dump.
The argument for compute displace/atom refresh is the ID of an
"atom-style variable"_variable.html which calculates a Boolean value
(0 or 1) based on the same criterion used by dump_modify thresh. What
this compute will then do is evaluate the atom-style variable and for
each atom that returns 1 (true), it will update the original
(reference) coordinates of the atom, which are stored by this compute.
The {refresh} argument for this compute is the ID of an "atom-style
variable"_variable.html which calculates a Boolean value (0 or 1)
based on the same criterion used by dump_modify thresh. This compute
evaluates the atom-style variable. For each atom that returns 1
(true), the original (reference) coordinates of the atom (stored by
this compute) are updated.
The effect of these commands is that a particular atom will only be
output in the dump file on shapshots following a diffusive hop. It
will not be output again until it performs another hop.
output in the dump file on the snapshot after it makes a diffusive
hop. It will not be output again until it makes another hop.
Note that the first snapshot (geneated by the first run after the dump
and this compute are defined), no atoms will be typically be output.
That is because the initial displacement for all atoms is 0.0. If an
initial dump snapshot is desired, containing the initial reference
position of all atoms, a command like this can be used before the
first simulation is performed
write_dump all custom tmp.dump.initial id type x y z :pre
Note that in the first snapshot of a subsequent run, no atoms will be
typically be output. That is because the initial displacement for all
atoms is 0.0. If an initial dump snapshot is desired, containing the
initial reference positions of all atoms, one way to do this is
illustrated above. An initial write_dump command can be used before
the first run. It will contain the positions of all the atoms,
Options in the "dump_modify"_dump_modify.html command above will
append new output to that same file and delay the output until a later
timestep. The {delay} setting avoids a second time = 0 snapshot which
would be empty.
:line

View File

@ -264,11 +264,15 @@ spacings.
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. When a simulation is performed on
different numbers of processors, there is no guarantee a particular
created atom will be assigned the same ID. If molecules are being
created, molecule IDs are assigned to created molecules in a similar
fashion.
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

View File

@ -80,7 +80,15 @@ deleted, then atom IDs are re-assigned so that they run from 1 to the
number of atoms in the system. Note that this is not done for
molecular systems (see the "atom_style"_atom_style.html command),
regardless of the {compress} setting, since it would foul up the bond
connectivity that has already been assigned.
connectivity that has already been assigned. However, the
"reset_ids"_reset_ids.html command can be used after this command to
accomplish the same thing.
Note that the re-assignement of IDs is not really a compression, where
gaps in atom IDs are removed by decrementing atom IDs that are larger.
Instead the IDs for all atoms are erased, and new IDs are assigned so
that the atoms owned by individual processors have consecutive IDs, as
the "create_atoms"_create_atoms.html command explains.
A molecular system with fixed bonds, angles, dihedrals, or improper
interactions, is one where the topology of the interactions is
@ -137,7 +145,7 @@ using molecule template files via the "molecule"_molecule.html and
[Related commands:]
"create_atoms"_create_atoms.html
"create_atoms"_create_atoms.html, "reset_ids"_reset_ids.html
[Default:]

View File

@ -121,8 +121,9 @@ See the "dump_modify every"_dump_modify.html command.
Only information for atoms in the specified group is dumped. The
"dump_modify thresh and region and refresh"_dump_modify.html commands
can also alter what atoms are included. Not all styles support all
these options; see details below.
can also alter what atoms are included. Not all styles support
these options; see details on the "dump_modify"_dump_modify.html doc
page.
As described below, the filename determines the kind of output (text
or binary or gzipped, one big file or one per timestep, one big file

View File

@ -15,11 +15,13 @@ dump_modify dump-ID keyword values ... :pre
dump-ID = ID of dump to modify :ulb,l
one or more keyword/value pairs may be appended :l
these keywords apply to various dump styles :l
keyword = {append} or {at} or {buffer} or {element} or {every} or {fileper} or {first} or {flush} or {format} or {image} or {label} or {nfile} or {pad} or {precision} or {region} or {scale} or {sort} or {thresh} or {unwrap} :l
keyword = {append} or {at} or {buffer} or {delay} or {element} or {every} or {fileper} or {first} or {flush} or {format} or {image} or {label} or {nfile} or {pad} or {precision} or {region} or {scale} or {sort} or {thresh} or {unwrap} :l
{append} arg = {yes} or {no}
{at} arg = N
N = index of frame written upon first dump
{buffer} arg = {yes} or {no}
{delay} arg = Dstep
Dstep = delay output until this timestep
{element} args = E1 E2 ... EN, where N = # of atom types
E1,...,EN = element name, e.g. C or Fe or Ga
{every} arg = N
@ -175,6 +177,14 @@ extra buffering.
:line
The {delay} keyword applies to all dump styles. No snapshots will be
output until the specified {Dstep} timestep or later. Specifying
{Dstep} < 0 is the same as turning off the delay setting. This is a
way to turn off unwanted output early in a simulation, for example,
during an equilibration phase.
:line
The {element} keyword applies only to the dump {cfg}, {xyz}, and
{image} styles. It associates element names (e.g. H, C, Fe) with
LAMMPS atom types. See the list of element names at the bottom of
@ -460,8 +470,11 @@ to respond to the call from the dump command, and update the
appropriate reference positions. This is done be defining an
"atom-style variable"_variable.html, {check} in this example, which
calculates a Boolean value (0 or 1) for each atom, based on the same
criterion used by dump_modify thresh. See the "compute
displace/atom"_compute_displace_atom.html command for more details.
criterion used by dump_modify thresh.
See the "compute displace/atom"_compute_displace_atom.html command for
more details, including an example of how to produce output that
includes an initial snapshot with the reference position of all atoms.
Note that only computes with a {refresh} option will work with
dump_modify refresh. See individual compute doc pages for details.

View File

@ -89,6 +89,7 @@ read_restart.html
region.html
replicate.html
rerun.html
reset_ids.html
reset_timestep.html
restart.html
run.html

56
doc/src/reset_ids.txt Normal file
View File

@ -0,0 +1,56 @@
"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
reset_ids command :h3
[Syntax:]
reset_ids :pre
[Examples:]
reset_ids :pre
[Description:]
Reset atom IDs for the system, including all the global IDs stored
for bond, angle, dihedral, improper topology data. This will
create a set of IDs that are numbered contiguously from 1 to N
for a N atoms system.
This can be useful to do after perfoming a "delete_atoms" command for
a molecular system. The delete_atoms compress yes option will not
perform this operation due to the existence of bond topology. It can
also be useful to do after any simulation which has lost atoms,
e.g. due to atoms moving outside a simulation box with fixed
boundaries (see the "boundary command"), or due to evaporation (see
the "fix evaporate" command).
Note that the resetting of IDs is not really a compression, where gaps
in atom IDs are removed by decrementing atom IDs that are larger.
Instead the IDs for all atoms are erased, and new IDs are assigned so
that the atoms owned by an individual processor have consecutive IDs,
as the "create_atoms"_create_atoms.html command explains.
NOTE: If this command is used before a "pair style"_pair_style.html is
defined, an error about bond topology atom IDs not being found may
result. This is because the cutoff distance for ghost atom
communication was not sufficient to find atoms in bonds, angles, etc
that are owned by other processors. The "comm_modify
cutoff"_comm_modify.html command can be used to correct this issue.
Or you can define a pair style before using this command. If you do
the former, you should unset the comm_modify cutoff after using
reset_ids so that subsequent communication is not inefficient.
[Restrictions:] none
[Related commands:]
"delete_atoms"_delete_atoms.html
[Default:] none

View File

@ -187,7 +187,8 @@ class Atom : protected Pointers {
int nextra_store;
int map_style; // style of atom map: 0=none, 1=array, 2=hash
int map_user; // user selected style = same 0,1,2
int map_user; // user requested map style:
// 0 = no request, 1=array, 2=hash, 3=yes
tagint map_tag_max; // max atom ID that map() is setup for
// spatial sorting of atoms

View File

@ -22,13 +22,13 @@
#include "force.h"
#include "group.h"
#include "region.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include <map>
@ -68,7 +68,8 @@ void DeleteAtoms::command(int narg, char **arg)
if (allflag) {
int igroup = group->find("all");
if ((igroup >= 0) && modify->check_rigid_group_overlap(group->bitmask[igroup]))
if ((igroup >= 0) &&
modify->check_rigid_group_overlap(group->bitmask[igroup]))
error->warning(FLERR,"Attempting to delete atoms in rigid bodies");
} else {
if (modify->check_rigid_list_overlap(dlist))
@ -105,7 +106,7 @@ void DeleteAtoms::command(int narg, char **arg)
memory->destroy(dlist);
}
// if non-molecular system and compress flag set,
// if non-molecular system and compress flag set:
// reset atom tags to be contiguous
// set all atom IDs to 0, call tag_extend()

View File

@ -89,6 +89,7 @@ Dump::Dump(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
buffer_flag = 0;
padflag = 0;
pbcflag = 0;
delay_flag = 0;
maxbuf = maxids = maxsort = maxproc = 0;
buf = bufsort = NULL;
@ -320,6 +321,10 @@ void Dump::write()
imageint *imagehold;
double **xhold,**vhold;
// if timestep < delaystep, just return
if (delay_flag && update->ntimestep < delaystep) return;
// if file per timestep, open new file
if (multifile) openfile();
@ -897,6 +902,13 @@ void Dump::modify_params(int narg, char **arg)
error->all(FLERR,"Dump_modify buffer yes not allowed for this style");
iarg += 2;
} else if (strcmp(arg[iarg],"delay") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
delaystep = force->bnumeric(FLERR,arg[iarg+1]);
if (delaystep >= 0) delay_flag = 1;
else delay_flag = 0;
iarg += 2;
} else if (strcmp(arg[iarg],"every") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump_modify command");
int idump;

View File

@ -77,6 +77,8 @@ class Dump : protected Pointers {
int sortcol; // 0 to sort on ID, 1-N on columns
int sortcolm1; // sortcol - 1
int sortorder; // ASCEND or DESCEND
int delay_flag; // 1 if delay output until delaystep
bigint delaystep;
int refreshflag; // 1 if dump_modify refresh specified
char *refresh; // compute ID to invoke refresh() on

View File

@ -577,13 +577,14 @@ int ReadDump::fields_and_keywords(int narg, char **arg)
fieldlabel = new char*[narg+2];
// add id and type fields as needed
// scan ahead to see if "add yes" keyword/value is used
// scan ahead to see if "add yes/keep" keyword/value is used
// requires extra "type" field from from dump file
int iarg;
for (iarg = 0; iarg < narg; iarg++)
if (strcmp(arg[iarg],"add") == 0)
if (iarg < narg-1 && strcmp(arg[iarg+1],"yes") == 0) break;
if (iarg < narg-1 && (strcmp(arg[iarg+1],"yes") == 0 ||
strcmp(arg[iarg+1],"keep") == 0)) break;
nfield = 0;
fieldtype[nfield++] = ID;

246
src/reset_ids.cpp Normal file
View File

@ -0,0 +1,246 @@
/* ----------------------------------------------------------------------
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 "reset_ids.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ResetIDs::ResetIDs(LAMMPS *lmp) : Pointers(lmp) {}
/* ---------------------------------------------------------------------- */
void ResetIDs::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Reset_ids command before simulation box is defined");
if (narg != 0) error->all(FLERR,"Illegal reset_ids command");
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use reset_ids unless atoms have IDs");
// NOTE: check if any fixes exist which store atom IDs?
// if so, this operation will mess up the fix
if (comm->me == 0) {
if (screen) fprintf(screen,"Resetting atom IDs ...\n");
if (logfile) fprintf(logfile,"Resetting atom IDs ...\n");
}
// create an atom map if one doesn't exist already
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// initialize system since comm->borders() will be invoked
lmp->init();
// setup domain, communication
// acquire ghosts - is that really necessary?
// exchange will clear map, borders will reset
// this is the map needed to lookup current global IDs for bond topology
if (domain->triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
comm->exchange();
comm->borders();
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
// oldIDs = copy of current owned IDs
tagint *tag = atom->tag;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
tagint *oldIDs;
memory->create(oldIDs,nlocal,"reset_ids:oldIDs");
for (int i = 0; i < nlocal; i++) {
oldIDs[i] = tag[i];
tag[i] = 0;
}
// assign new contigous IDs to owned atoms via tag_extend()
atom->tag_extend();
// newIDs = copy of new IDs
// restore old IDs, consistent with existing atom map
// forward_comm_array acquires new IDs for ghost atoms
double **newIDs;
memory->create(newIDs,nall,1,"reset_ids:newIDs");
for (int i = 0; i < nlocal; i++) {
newIDs[i][0] = tag[i];
tag[i] = oldIDs[i];
}
comm->forward_comm_array(1,newIDs);
// loop over bonds, angles, etc and reset IDs in stored topology arrays
// only necessary for molecular = 1, not molecular = 2
// badcount = atom IDs that could not be found
int badcount = 0;
if (atom->molecular == 1) {
int j,m;
tagint oldID;
if (atom->avec->bonds_allow) {
int *num_bond = atom->num_bond;
tagint **bond_atom = atom->bond_atom;
for (int i = 0; i < nlocal; i++) {
for (j = 0; j < num_bond[i]; j++) {
oldID = bond_atom[i][j];
m = atom->map(oldID);
if (m >= 0) bond_atom[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
}
}
}
if (atom->avec->angles_allow) {
int *num_angle = atom->num_angle;
tagint **angle_atom1 = atom->angle_atom1;
tagint **angle_atom2 = atom->angle_atom2;
tagint **angle_atom3 = atom->angle_atom3;
for (int i = 0; i < nlocal; i++) {
for (j = 0; j < num_angle[i]; j++) {
oldID = angle_atom1[i][j];
m = atom->map(oldID);
if (m >= 0) angle_atom1[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = angle_atom2[i][j];
m = atom->map(oldID);
if (m >= 0) angle_atom2[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = angle_atom3[i][j];
m = atom->map(oldID);
if (m >= 0) angle_atom3[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
}
}
}
if (atom->avec->dihedrals_allow) {
int *num_dihedral = atom->num_dihedral;
tagint **dihedral_atom1 = atom->dihedral_atom1;
tagint **dihedral_atom2 = atom->dihedral_atom2;
tagint **dihedral_atom3 = atom->dihedral_atom3;
tagint **dihedral_atom4 = atom->dihedral_atom4;
for (int i = 0; i < nlocal; i++) {
for (j = 0; j < num_dihedral[i]; j++) {
oldID = dihedral_atom1[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom1[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = dihedral_atom2[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom2[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = dihedral_atom3[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom3[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = dihedral_atom4[i][j];
m = atom->map(oldID);
if (m >= 0) dihedral_atom4[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
}
}
}
if (atom->avec->impropers_allow) {
int *num_improper = atom->num_improper;
tagint **improper_atom1 = atom->improper_atom1;
tagint **improper_atom2 = atom->improper_atom2;
tagint **improper_atom3 = atom->improper_atom3;
tagint **improper_atom4 = atom->improper_atom4;
for (int i = 0; i < nlocal; i++) {
for (j = 0; j < num_improper[i]; j++) {
oldID = improper_atom1[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom1[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = improper_atom2[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom2[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = improper_atom3[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom3[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
oldID = improper_atom4[i][j];
m = atom->map(oldID);
if (m >= 0) improper_atom4[i][j] = static_cast<tagint> (newIDs[m][0]);
else badcount++;
}
}
}
}
// error check
int all;
MPI_Allreduce(&badcount,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,"Reset_ids missing %d bond topology atom IDs - "
"use comm_modify cutoff",all);
error->all(FLERR,str);
}
// reset IDs and atom map for owned atoms
atom->map_clear();
atom->nghost = 0;
for (int i = 0; i < nlocal; i++) tag[i] = static_cast<tagint> (newIDs[i][0]);
atom->map_init();
atom->map_set();
// delete temporary atom map
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
}
// clean up
memory->destroy(oldIDs);
memory->destroy(newIDs);
}

42
src/reset_ids.h Normal file
View File

@ -0,0 +1,42 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifdef COMMAND_CLASS
CommandStyle(reset_ids,ResetIDs)
#else
#ifndef LMP_RESET_IDS_H
#define LMP_RESET_IDS_H
#include "pointers.h"
namespace LAMMPS_NS {
class ResetIDs : protected Pointers {
public:
ResetIDs(class LAMMPS *);
void command(int, char **);
private:
};
}
#endif
#endif
/* ERROR/WARNING messages:
*/