diff --git a/src/atom.cpp b/src/atom.cpp
index ac88c4308..e96e5e205 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1,2166 +1,2185 @@
/* ----------------------------------------------------------------------
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 <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include "atom.h"
#include "style_atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "neighbor.h"
#include "force.h"
#include "modify.h"
#include "fix.h"
#include "output.h"
#include "thermo.h"
#include "update.h"
#include "domain.h"
#include "group.h"
#include "molecule.h"
#include "accelerator_cuda.h"
#include "atom_masks.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 1
#define DELTA_MEMSTR 1024
#define EPSILON 1.0e-6
#define CUDA_CHUNK 3000
enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
/* ---------------------------------------------------------------------- */
Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
{
natoms = 0;
nlocal = nghost = nmax = 0;
ntypes = 0;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
nbonds = nangles = ndihedrals = nimpropers = 0;
firstgroupname = NULL;
sortfreq = 1000;
nextsort = 0;
userbinsize = 0.0;
maxbin = maxnext = 0;
binhead = NULL;
next = permute = NULL;
// initialize atom arrays
// customize by adding new array
tag = NULL;
type = mask = NULL;
image = NULL;
x = v = f = NULL;
molecule = NULL;
molindex = molatom = NULL;
q = NULL;
mu = NULL;
omega = angmom = torque = NULL;
radius = rmass = NULL;
ellipsoid = line = tri = body = NULL;
vfrac = s0 = NULL;
x0 = NULL;
spin = NULL;
eradius = ervel = erforce = NULL;
cs = csforce = vforce = ervelforce = NULL;
etag = NULL;
rho = drho = e = de = cv = NULL;
vest = NULL;
+ // USER-DPD
+
+ uCond = uMech = uChem = uCG = uCGnew = NULL;
+ duCond = duMech = duChem = NULL;
+ dpdTheta = NULL;
+
// USER-SMD
contact_radius = NULL;
smd_data_9 = NULL;
smd_stress = NULL;
eff_plastic_strain = NULL;
eff_plastic_strain_rate = NULL;
damage = NULL;
// molecular info
bond_per_atom = extra_bond_per_atom = 0;
num_bond = NULL;
bond_type = NULL;
bond_atom = NULL;
angle_per_atom = extra_angle_per_atom = 0;
num_angle = NULL;
angle_type = NULL;
angle_atom1 = angle_atom2 = angle_atom3 = NULL;
dihedral_per_atom = extra_dihedral_per_atom = 0;
num_dihedral = NULL;
dihedral_type = NULL;
dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = NULL;
improper_per_atom = extra_improper_per_atom = 0;
num_improper = NULL;
improper_type = NULL;
improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = NULL;
maxspecial = 1;
nspecial = NULL;
special = NULL;
// user-defined molecules
nmolecule = 0;
molecules = NULL;
// custom atom arrays
nivector = ndvector = 0;
ivector = NULL;
dvector = NULL;
iname = dname = NULL;
// initialize atom style and array existence flags
// customize by adding new flag
sphere_flag = peri_flag = electron_flag = 0;
wavepacket_flag = sph_flag = 0;
molecule_flag = 0;
q_flag = mu_flag = 0;
omega_flag = torque_flag = angmom_flag = 0;
radius_flag = rmass_flag = 0;
ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
vfrac_flag = 0;
spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0;
cs_flag = csforce_flag = vforce_flag = etag_flag = 0;
rho_flag = e_flag = cv_flag = vest_flag = 0;
+ dpd_flag = 0;
// USER-SMD
smd_flag = 0;
contact_radius_flag = 0;
smd_data_9_flag = 0;
smd_stress_flag = 0;
x0_flag = 0;
eff_plastic_strain_flag = 0;
eff_plastic_strain_rate_flag = 0;
damage_flag = 0;
// Peridynamic scale factor
pdscale = 1.0;
// ntype-length arrays
mass = NULL;
mass_setflag = NULL;
// callback lists & extra restart info
nextra_grow = nextra_restart = nextra_border = 0;
extra_grow = extra_restart = extra_border = NULL;
nextra_grow_max = nextra_restart_max = nextra_border_max = 0;
nextra_store = 0;
extra = NULL;
// default atom ID and mapping values
tag_enable = 1;
map_style = map_user = 0;
map_tag_max = -1;
map_maxarray = map_nhash = -1;
max_same = 0;
sametag = NULL;
map_array = NULL;
map_bucket = NULL;
map_hash = NULL;
atom_style = NULL;
avec = NULL;
datamask = ALL_MASK;
datamask_ext = ALL_MASK;
}
/* ---------------------------------------------------------------------- */
Atom::~Atom()
{
delete [] atom_style;
delete avec;
delete [] firstgroupname;
memory->destroy(binhead);
memory->destroy(next);
memory->destroy(permute);
// delete atom arrays
// customize by adding new array
memory->destroy(tag);
memory->destroy(type);
memory->destroy(mask);
memory->destroy(image);
memory->destroy(x);
memory->destroy(v);
memory->destroy(f);
memory->destroy(molecule);
memory->destroy(molindex);
memory->destroy(molatom);
memory->destroy(q);
memory->destroy(mu);
memory->destroy(omega);
memory->destroy(angmom);
memory->destroy(torque);
memory->destroy(radius);
memory->destroy(rmass);
memory->destroy(ellipsoid);
memory->destroy(line);
memory->destroy(tri);
memory->destroy(body);
memory->destroy(vfrac);
memory->destroy(s0);
memory->destroy(x0);
memory->destroy(spin);
memory->destroy(eradius);
memory->destroy(ervel);
memory->destroy(erforce);
memory->destroy(ervelforce);
memory->destroy(cs);
memory->destroy(csforce);
memory->destroy(vforce);
memory->destroy(etag);
memory->destroy(rho);
memory->destroy(drho);
memory->destroy(e);
memory->destroy(de);
memory->destroy(cv);
memory->destroy(vest);
memory->destroy(contact_radius);
memory->destroy(smd_data_9);
memory->destroy(smd_stress);
memory->destroy(eff_plastic_strain);
memory->destroy(eff_plastic_strain_rate);
memory->destroy(damage);
+ memory->destroy(dpdTheta);
+ memory->destroy(uCond);
+ memory->destroy(uMech);
+ memory->destroy(uChem);
+ memory->destroy(uCG);
+ memory->destroy(uCGnew);
+ memory->destroy(duCond);
+ memory->destroy(duMech);
+ memory->destroy(duChem);
+
memory->destroy(nspecial);
memory->destroy(special);
memory->destroy(num_bond);
memory->destroy(bond_type);
memory->destroy(bond_atom);
memory->destroy(num_angle);
memory->destroy(angle_type);
memory->destroy(angle_atom1);
memory->destroy(angle_atom2);
memory->destroy(angle_atom3);
memory->destroy(num_dihedral);
memory->destroy(dihedral_type);
memory->destroy(dihedral_atom1);
memory->destroy(dihedral_atom2);
memory->destroy(dihedral_atom3);
memory->destroy(dihedral_atom4);
memory->destroy(num_improper);
memory->destroy(improper_type);
memory->destroy(improper_atom1);
memory->destroy(improper_atom2);
memory->destroy(improper_atom3);
memory->destroy(improper_atom4);
// delete custom atom arrays
for (int i = 0; i < nivector; i++) {
delete [] iname[i];
memory->destroy(ivector[i]);
}
for (int i = 0; i < ndvector; i++) {
delete [] dname[i];
memory->destroy(dvector[i]);
}
memory->sfree(iname);
memory->sfree(dname);
memory->sfree(ivector);
memory->sfree(dvector);
// delete user-defined molecules
for (int i = 0; i < nmolecule; i++) delete molecules[i];
memory->sfree(molecules);
// delete per-type arrays
delete [] mass;
delete [] mass_setflag;
// delete extra arrays
memory->destroy(extra_grow);
memory->destroy(extra_restart);
memory->destroy(extra_border);
memory->destroy(extra);
// delete mapping data structures
map_delete();
}
/* ----------------------------------------------------------------------
copy modify settings from old Atom class to current Atom class
------------------------------------------------------------------------- */
void Atom::settings(Atom *old)
{
tag_enable = old->tag_enable;
map_user = old->map_user;
map_style = old->map_style;
sortfreq = old->sortfreq;
userbinsize = old->userbinsize;
if (old->firstgroupname) {
int n = strlen(old->firstgroupname) + 1;
firstgroupname = new char[n];
strcpy(firstgroupname,old->firstgroupname);
}
}
/* ----------------------------------------------------------------------
create an AtomVec style
called from lammps.cpp, input script, restart file, replicate
------------------------------------------------------------------------- */
void Atom::create_avec(const char *style, int narg, char **arg, int trysuffix)
{
delete [] atom_style;
if (avec) delete avec;
// unset atom style and array existence flags
// may have been set by old avec
// customize by adding new flag
sphere_flag = peri_flag = electron_flag = 0;
wavepacket_flag = sph_flag = 0;
molecule_flag = 0;
q_flag = mu_flag = 0;
omega_flag = torque_flag = angmom_flag = 0;
radius_flag = rmass_flag = 0;
ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
vfrac_flag = 0;
spin_flag = eradius_flag = ervel_flag = erforce_flag = ervelforce_flag = 0;
cs_flag = csforce_flag = vforce_flag = etag_flag = 0;
rho_flag = e_flag = cv_flag = vest_flag = 0;
// create instance of AtomVec
// use grow() to initialize atom-based arrays to length 1
// so that x[0][0] can always be referenced even if proc has no atoms
int sflag;
avec = new_avec(style,trysuffix,sflag);
avec->store_args(narg,arg);
avec->process_args(narg,arg);
avec->grow(1);
if (sflag) {
char estyle[256];
if (sflag == 1) sprintf(estyle,"%s/%s",style,lmp->suffix);
else sprintf(estyle,"%s/%s",style,lmp->suffix2);
int n = strlen(estyle) + 1;
atom_style = new char[n];
strcpy(atom_style,estyle);
} else {
int n = strlen(style) + 1;
atom_style = new char[n];
strcpy(atom_style,style);
}
// if molecular system:
// atom IDs must be defined
// force atom map to be created
// map style may be reset by map_init() and its call to map_style_set()
molecular = avec->molecular;
if (molecular && tag_enable == 0)
error->all(FLERR,"Atom IDs must be used for molecular systems");
if (molecular) map_style = 1;
}
/* ----------------------------------------------------------------------
generate an AtomVec class, first with suffix appended
------------------------------------------------------------------------- */
AtomVec *Atom::new_avec(const char *style, int trysuffix, int &sflag)
{
if (trysuffix && lmp->suffix_enable) {
if (lmp->suffix) {
sflag = 1;
char estyle[256];
sprintf(estyle,"%s/%s",style,lmp->suffix);
if (0) return NULL;
#define ATOM_CLASS
#define AtomStyle(key,Class) \
else if (strcmp(estyle,#key) == 0) return new Class(lmp);
#include "style_atom.h"
#undef AtomStyle
#undef ATOM_CLASS
}
if (lmp->suffix2) {
sflag = 2;
char estyle[256];
sprintf(estyle,"%s/%s",style,lmp->suffix2);
if (0) return NULL;
#define ATOM_CLASS
#define AtomStyle(key,Class) \
else if (strcmp(estyle,#key) == 0) return new Class(lmp);
#include "style_atom.h"
#undef AtomStyle
#undef ATOM_CLASS
}
}
sflag = 0;
if (0) return NULL;
#define ATOM_CLASS
#define AtomStyle(key,Class) \
else if (strcmp(style,#key) == 0) return new Class(lmp);
#include "style_atom.h"
#undef ATOM_CLASS
else error->all(FLERR,"Unknown atom style");
return NULL;
}
/* ---------------------------------------------------------------------- */
void Atom::init()
{
// delete extra array since it doesn't persist past first run
if (nextra_store) {
memory->destroy(extra);
extra = NULL;
nextra_store = 0;
}
// check arrays that are atom type in length
check_mass();
// setup of firstgroup
if (firstgroupname) {
firstgroup = group->find(firstgroupname);
if (firstgroup < 0)
error->all(FLERR,"Could not find atom_modify first group ID");
} else firstgroup = -1;
// init AtomVec
avec->init();
}
/* ---------------------------------------------------------------------- */
void Atom::setup()
{
// setup bins for sorting
// cannot do this in init() because uses neighbor cutoff
if (sortfreq > 0) setup_sort_bins();
}
/* ----------------------------------------------------------------------
return ptr to AtomVec class if matches style or to matching hybrid sub-class
return NULL if no match
------------------------------------------------------------------------- */
AtomVec *Atom::style_match(const char *style)
{
if (strcmp(atom_style,style) == 0) return avec;
else if (strcmp(atom_style,"hybrid") == 0) {
AtomVecHybrid *avec_hybrid = (AtomVecHybrid *) avec;
for (int i = 0; i < avec_hybrid->nstyles; i++)
if (strcmp(avec_hybrid->keywords[i],style) == 0)
return avec_hybrid->styles[i];
}
return NULL;
}
/* ----------------------------------------------------------------------
modify parameters of the atom style
some options can only be invoked before simulation box is defined
first and sort options cannot be used together
------------------------------------------------------------------------- */
void Atom::modify_params(int narg, char **arg)
{
if (narg == 0) error->all(FLERR,"Illegal atom_modify command");
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"id") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
if (domain->box_exist)
error->all(FLERR,
"Atom_modify id command after simulation box is defined");
if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1;
else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 0;
else error->all(FLERR,"Illegal atom_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"map") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
if (domain->box_exist)
error->all(FLERR,
"Atom_modify map command after simulation box is defined");
if (strcmp(arg[iarg+1],"array") == 0) map_user = 1;
else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2;
else error->all(FLERR,"Illegal atom_modify command");
map_style = map_user;
iarg += 2;
} else if (strcmp(arg[iarg],"first") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
if (strcmp(arg[iarg+1],"all") == 0) {
delete [] firstgroupname;
firstgroupname = NULL;
} else {
int n = strlen(arg[iarg+1]) + 1;
firstgroupname = new char[n];
strcpy(firstgroupname,arg[iarg+1]);
sortfreq = 0;
}
iarg += 2;
} else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command");
sortfreq = force->inumeric(FLERR,arg[iarg+1]);
userbinsize = force->numeric(FLERR,arg[iarg+2]);
if (sortfreq < 0 || userbinsize < 0.0)
error->all(FLERR,"Illegal atom_modify command");
if (sortfreq >= 0 && firstgroupname)
error->all(FLERR,"Atom_modify sort and first options "
"cannot be used together");
iarg += 3;
} else error->all(FLERR,"Illegal atom_modify command");
}
}
/* ----------------------------------------------------------------------
check that atom IDs are valid
error if any atom ID < 0 or atom ID = MAXTAGINT
if any atom ID > 0, error if any atom ID == 0
if any atom ID > 0, error if tag_enable = 0
if all atom IDs = 0, tag_enable must be 0
if max atom IDs < natoms, must be duplicates
OK if max atom IDs > natoms
NOTE: not fully checking that atom IDs are unique
------------------------------------------------------------------------- */
void Atom::tag_check()
{
tagint min = MAXTAGINT;
tagint max = 0;
for (int i = 0; i < nlocal; i++) {
min = MIN(min,tag[i]);
max = MAX(max,tag[i]);
}
tagint minall,maxall;
MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world);
MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
if (minall < 0) error->all(FLERR,"One or more Atom IDs is negative");
if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs is too big");
if (maxall > 0 && minall == 0)
error->all(FLERR,"One or more atom IDs is zero");
if (maxall > 0 && tag_enable == 0)
error->all(FLERR,"Non-zero atom IDs with atom_modify id = no");
if (maxall == 0 && natoms && tag_enable)
error->all(FLERR,"All atom IDs = 0 but atom_modify id = yes");
if (tag_enable && maxall < natoms)
error->all(FLERR,"Duplicate atom IDs exist");
}
/* ----------------------------------------------------------------------
add unique tags to any atoms with tag = 0
new tags are grouped by proc and start after max current tag
called after creating new atoms
error if new tags will exceed MAXTAGINT
------------------------------------------------------------------------- */
void Atom::tag_extend()
{
// maxtag_all = max tag for all atoms
tagint maxtag = 0;
for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
tagint maxtag_all;
MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
// DEBUG: useful for generating 64-bit IDs even for small systems
// use only when LAMMPS is compiled with BIGBIG
//maxtag_all += 1000000000000;
// notag = # of atoms I own with no tag (tag = 0)
// notag_sum = # of total atoms on procs <= me with no tag
bigint notag = 0;
for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++;
bigint notag_total;
MPI_Allreduce(¬ag,¬ag_total,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (notag_total >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
bigint notag_sum;
MPI_Scan(¬ag,¬ag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
// itag = 1st new tag that my untagged atoms should use
tagint itag = maxtag_all + notag_sum - notag + 1;
for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++;
}
/* ----------------------------------------------------------------------
check that atom IDs span range from 1 to Natoms inclusive
return 0 if mintag != 1 or maxtag != Natoms
return 1 if OK
doesn't actually check if all tag values are used
------------------------------------------------------------------------- */
int Atom::tag_consecutive()
{
tagint idmin = MAXTAGINT;
tagint idmax = 0;
for (int i = 0; i < nlocal; i++) {
idmin = MIN(idmin,tag[i]);
idmax = MAX(idmax,tag[i]);
}
tagint idminall,idmaxall;
MPI_Allreduce(&idmin,&idminall,1,MPI_LMP_TAGINT,MPI_MIN,world);
MPI_Allreduce(&idmax,&idmaxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
if (idminall != 1 || idmaxall != natoms) return 0;
return 1;
}
/* ----------------------------------------------------------------------
count and return words in a single line
make copy of line before using strtok so as not to change line
trim anything from '#' onward
------------------------------------------------------------------------- */
int Atom::count_words(const char *line)
{
int n = strlen(line) + 1;
char *copy;
memory->create(copy,n,"atom:copy");
strcpy(copy,line);
char *ptr;
if ((ptr = strchr(copy,'#'))) *ptr = '\0';
if (strtok(copy," \t\n\r\f") == NULL) {
memory->destroy(copy);
return 0;
}
n = 1;
while (strtok(NULL," \t\n\r\f")) n++;
memory->destroy(copy);
return n;
}
/* ----------------------------------------------------------------------
count and return words in a single line using provided copy buf
make copy of line before using strtok so as not to change line
trim anything from '#' onward
------------------------------------------------------------------------- */
int Atom::count_words(const char *line, char *copy)
{
strcpy(copy,line);
char *ptr;
if ((ptr = strchr(copy,'#'))) *ptr = '\0';
if (strtok(copy," \t\n\r\f") == NULL) {
memory->destroy(copy);
return 0;
}
int n = 1;
while (strtok(NULL," \t\n\r\f")) n++;
return n;
}
/* ----------------------------------------------------------------------
deallocate molecular topology arrays
done before realloc with (possibly) new 2nd dimension set to
correctly initialized per-atom values, e.g. bond_per_atom
needs to be called whenever 2nd dimensions are changed
and these arrays are already pre-allocated,
e.g. due to grow(1) in create_avec()
------------------------------------------------------------------------- */
void Atom::deallocate_topology()
{
memory->destroy(atom->bond_type);
memory->destroy(atom->bond_atom);
atom->bond_type = NULL;
atom->bond_atom = NULL;
memory->destroy(atom->angle_type);
memory->destroy(atom->angle_atom1);
memory->destroy(atom->angle_atom2);
memory->destroy(atom->angle_atom3);
atom->angle_type = NULL;
atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL;
memory->destroy(atom->dihedral_type);
memory->destroy(atom->dihedral_atom1);
memory->destroy(atom->dihedral_atom2);
memory->destroy(atom->dihedral_atom3);
memory->destroy(atom->dihedral_atom4);
atom->dihedral_type = NULL;
atom->dihedral_atom1 = atom->dihedral_atom2 =
atom->dihedral_atom3 = atom->dihedral_atom4 = NULL;
memory->destroy(atom->improper_type);
memory->destroy(atom->improper_atom1);
memory->destroy(atom->improper_atom2);
memory->destroy(atom->improper_atom3);
memory->destroy(atom->improper_atom4);
atom->improper_type = NULL;
atom->improper_atom1 = atom->improper_atom2 =
atom->improper_atom3 = atom->improper_atom4 = NULL;
}
/* ----------------------------------------------------------------------
unpack N lines from Atom section of data file
call style-specific routine to parse line
------------------------------------------------------------------------- */
void Atom::data_atoms(int n, char *buf, tagint id_offset, int type_offset,
int shiftflag, double *shift)
{
int m,xptr,iptr;
imageint imagedata;
double xdata[3],lamda[3];
double *coord;
char *next;
next = strchr(buf,'\n');
*next = '\0';
int nwords = count_words(buf);
*next = '\n';
if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3)
error->all(FLERR,"Incorrect atom format in data file");
char **values = new char*[nwords];
// set bounds for my proc
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
// insures all data atoms will be owned even with round-off
int triclinic = domain->triclinic;
double epsilon[3];
if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
else {
epsilon[0] = domain->prd[0] * EPSILON;
epsilon[1] = domain->prd[1] * EPSILON;
epsilon[2] = domain->prd[2] * EPSILON;
}
double sublo[3],subhi[3];
if (triclinic == 0) {
sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
} else {
sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
}
if (comm->layout != LAYOUT_TILED) {
if (domain->xperiodic) {
if (comm->myloc[0] == 0) sublo[0] -= epsilon[0];
if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0];
}
if (domain->yperiodic) {
if (comm->myloc[1] == 0) sublo[1] -= epsilon[1];
if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1];
}
if (domain->zperiodic) {
if (comm->myloc[2] == 0) sublo[2] -= epsilon[2];
if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2];
}
} else {
if (domain->xperiodic) {
if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0];
if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0];
}
if (domain->yperiodic) {
if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1];
if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1];
}
if (domain->zperiodic) {
if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2];
if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2];
}
}
// xptr = which word in line starts xyz coords
// iptr = which word in line starts ix,iy,iz image flags
xptr = avec->xcol_data - 1;
int imageflag = 0;
if (nwords > avec->size_data_atom) imageflag = 1;
if (imageflag) iptr = nwords - 3;
// loop over lines of atom data
// tokenize the line into values
// extract xyz coords and image flags
// remap atom into simulation box
// if atom is in my sub-domain, unpack its values
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
if (values[0] == NULL)
error->all(FLERR,"Incorrect atom format in data file");
for (m = 1; m < nwords; m++) {
values[m] = strtok(NULL," \t\n\r\f");
if (values[m] == NULL)
error->all(FLERR,"Incorrect atom format in data file");
}
if (imageflag)
imagedata = ((imageint) (atoi(values[iptr]) + IMGMAX) & IMGMASK) |
(((imageint) (atoi(values[iptr+1]) + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS);
else imagedata = ((imageint) IMGMAX << IMG2BITS) |
((imageint) IMGMAX << IMGBITS) | IMGMAX;
xdata[0] = atof(values[xptr]);
xdata[1] = atof(values[xptr+1]);
xdata[2] = atof(values[xptr+2]);
if (shiftflag) {
xdata[0] += shift[0];
xdata[1] += shift[1];
xdata[2] += shift[2];
}
domain->remap(xdata,imagedata);
if (triclinic) {
domain->x2lamda(xdata,lamda);
coord = lamda;
} else coord = xdata;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) {
avec->data_atom(xdata,imagedata,values);
if (id_offset) tag[nlocal-1] += id_offset;
if (type_offset) {
type[nlocal-1] += type_offset;
if (type[nlocal-1] > ntypes)
error->one(FLERR,"Invalid atom type in Atoms section of data file");
}
}
buf = next + 1;
}
delete [] values;
}
/* ----------------------------------------------------------------------
unpack N lines from Velocity section of data file
check that atom IDs are > 0 and <= map_tag_max
call style-specific routine to parse line
------------------------------------------------------------------------- */
void Atom::data_vels(int n, char *buf, tagint id_offset)
{
int j,m;
tagint tagdata;
char *next;
next = strchr(buf,'\n');
*next = '\0';
int nwords = count_words(buf);
*next = '\n';
if (nwords != avec->size_data_vel)
error->all(FLERR,"Incorrect velocity format in data file");
char **values = new char*[nwords];
// loop over lines of atom velocities
// tokenize the line into values
// if I own atom tag, unpack its values
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
for (j = 1; j < nwords; j++)
values[j] = strtok(NULL," \t\n\r\f");
tagdata = ATOTAGINT(values[0]) + id_offset;
if (tagdata <= 0 || tagdata > map_tag_max)
error->one(FLERR,"Invalid atom ID in Velocities section of data file");
if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]);
buf = next + 1;
}
delete [] values;
}
/* ----------------------------------------------------------------------
process N bonds read into buf from data files
if count is non-NULL, just count bonds per atom
else store them with atoms
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset,
int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2;
char *next;
int newton_bond = force->newton_bond;
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
*next = '\0';
sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
}
itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max)
error->one(FLERR,"Invalid atom ID in Bonds section of data file");
if (itype <= 0 || itype > nbondtypes)
error->one(FLERR,"Invalid bond type in Bonds section of data file");
if ((m = map(atom1)) >= 0) {
if (count) count[m]++;
else {
bond_type[m][num_bond[m]] = itype;
bond_atom[m][num_bond[m]] = atom2;
num_bond[m]++;
}
}
if (newton_bond == 0) {
if ((m = map(atom2)) >= 0) {
if (count) count[m]++;
else {
bond_type[m][num_bond[m]] = itype;
bond_atom[m][num_bond[m]] = atom1;
num_bond[m]++;
}
}
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
process N angles read into buf from data files
if count is non-NULL, just count angles per atom
else store them with atoms
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
void Atom::data_angles(int n, char *buf, int *count, tagint id_offset,
int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2,atom3;
char *next;
int newton_bond = force->newton_bond;
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
*next = '\0';
sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2,&atom3);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
atom3 += id_offset;
}
itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max ||
atom3 <= 0 || atom3 > map_tag_max)
error->one(FLERR,"Invalid atom ID in Angles section of data file");
if (itype <= 0 || itype > nangletypes)
error->one(FLERR,"Invalid angle type in Angles section of data file");
if ((m = map(atom2)) >= 0) {
if (count) count[m]++;
else {
angle_type[m][num_angle[m]] = itype;
angle_atom1[m][num_angle[m]] = atom1;
angle_atom2[m][num_angle[m]] = atom2;
angle_atom3[m][num_angle[m]] = atom3;
num_angle[m]++;
}
}
if (newton_bond == 0) {
if ((m = map(atom1)) >= 0) {
if (count) count[m]++;
else {
angle_type[m][num_angle[m]] = itype;
angle_atom1[m][num_angle[m]] = atom1;
angle_atom2[m][num_angle[m]] = atom2;
angle_atom3[m][num_angle[m]] = atom3;
num_angle[m]++;
}
}
if ((m = map(atom3)) >= 0) {
if (count) count[m]++;
else {
angle_type[m][num_angle[m]] = itype;
angle_atom1[m][num_angle[m]] = atom1;
angle_atom2[m][num_angle[m]] = atom2;
angle_atom3[m][num_angle[m]] = atom3;
num_angle[m]++;
}
}
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
process N dihedrals read into buf from data files
if count is non-NULL, just count diihedrals per atom
else store them with atoms
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset,
int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2,atom3,atom4;
char *next;
int newton_bond = force->newton_bond;
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
*next = '\0';
sscanf(buf,"%d %d "
TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2,&atom3,&atom4);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
atom3 += id_offset;
atom4 += id_offset;
}
itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max ||
atom3 <= 0 || atom3 > map_tag_max ||
atom4 <= 0 || atom4 > map_tag_max)
error->one(FLERR,"Invalid atom ID in Dihedrals section of data file");
if (itype <= 0 || itype > ndihedraltypes)
error->one(FLERR,
"Invalid dihedral type in Dihedrals section of data file");
if ((m = map(atom2)) >= 0) {
if (count) count[m]++;
else {
dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++;
}
}
if (newton_bond == 0) {
if ((m = map(atom1)) >= 0) {
if (count) count[m]++;
else {
dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++;
}
}
if ((m = map(atom3)) >= 0) {
if (count) count[m]++;
else {
dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++;
}
}
if ((m = map(atom4)) >= 0) {
if (count) count[m]++;
else {
dihedral_type[m][num_dihedral[m]] = itype;
dihedral_atom1[m][num_dihedral[m]] = atom1;
dihedral_atom2[m][num_dihedral[m]] = atom2;
dihedral_atom3[m][num_dihedral[m]] = atom3;
dihedral_atom4[m][num_dihedral[m]] = atom4;
num_dihedral[m]++;
}
}
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
process N impropers read into buf from data files
if count is non-NULL, just count impropers per atom
else store them with atoms
check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */
void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset,
int type_offset)
{
int m,tmp,itype;
tagint atom1,atom2,atom3,atom4;
char *next;
int newton_bond = force->newton_bond;
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
*next = '\0';
sscanf(buf,"%d %d "
TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
&tmp,&itype,&atom1,&atom2,&atom3,&atom4);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
atom3 += id_offset;
atom4 += id_offset;
}
itype += type_offset;
if (atom1 <= 0 || atom1 > map_tag_max ||
atom2 <= 0 || atom2 > map_tag_max ||
atom3 <= 0 || atom3 > map_tag_max ||
atom4 <= 0 || atom4 > map_tag_max)
error->one(FLERR,"Invalid atom ID in Impropers section of data file");
if (itype <= 0 || itype > nimpropertypes)
error->one(FLERR,
"Invalid improper type in Impropers section of data file");
if ((m = map(atom2)) >= 0) {
if (count) count[m]++;
else {
improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++;
}
}
if (newton_bond == 0) {
if ((m = map(atom1)) >= 0) {
if (count) count[m]++;
else {
improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++;
}
}
if ((m = map(atom3)) >= 0) {
if (count) count[m]++;
else {
improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++;
}
}
if ((m = map(atom4)) >= 0) {
if (count) count[m]++;
else {
improper_type[m][num_improper[m]] = itype;
improper_atom1[m][num_improper[m]] = atom1;
improper_atom2[m][num_improper[m]] = atom2;
improper_atom3[m][num_improper[m]] = atom3;
improper_atom4[m][num_improper[m]] = atom4;
num_improper[m]++;
}
}
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
unpack N lines from atom-style specific bonus section of data file
check that atom IDs are > 0 and <= map_tag_max
call style-specific routine to parse line
------------------------------------------------------------------------- */
void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset)
{
int j,m,tagdata;
char *next;
next = strchr(buf,'\n');
*next = '\0';
int nwords = count_words(buf);
*next = '\n';
if (nwords != avec_bonus->size_data_bonus)
error->all(FLERR,"Incorrect bonus data format in data file");
char **values = new char*[nwords];
// loop over lines of bonus atom data
// tokenize the line into values
// if I own atom tag, unpack its values
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
for (j = 1; j < nwords; j++)
values[j] = strtok(NULL," \t\n\r\f");
tagdata = ATOTAGINT(values[0]) + id_offset;
if (tagdata <= 0 || tagdata > map_tag_max)
error->one(FLERR,"Invalid atom ID in Bonus section of data file");
// ok to call child's data_atom_bonus() method thru parent avec_bonus,
// since data_bonus() was called with child ptr, and method is virtual
if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]);
buf = next + 1;
}
delete [] values;
}
/* ----------------------------------------------------------------------
unpack N bodies from Bodies section of data file
each body spans multiple lines
check that atom IDs are > 0 and <= map_tag_max
call style-specific routine to parse line
------------------------------------------------------------------------- */
void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body,
tagint id_offset)
{
int j,m,nvalues,tagdata,ninteger,ndouble;
int maxint = 0;
int maxdouble = 0;
int *ivalues = NULL;
double *dvalues = NULL;
// loop over lines of body data
// if I own atom tag, tokenize lines into ivalues/dvalues, call data_body()
// else skip values
for (int i = 0; i < n; i++) {
if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset;
else tagdata = ATOTAGINT(strtok(NULL," \t\n\r\f")) + id_offset;
if (tagdata <= 0 || tagdata > map_tag_max)
error->one(FLERR,"Invalid atom ID in Bodies section of data file");
ninteger = force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
ndouble = force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
if ((m = map(tagdata)) >= 0) {
if (ninteger > maxint) {
delete [] ivalues;
maxint = ninteger;
ivalues = new int[maxint];
}
if (ndouble > maxdouble) {
delete [] dvalues;
maxdouble = ndouble;
dvalues = new double[maxdouble];
}
for (j = 0; j < ninteger; j++)
ivalues[j] = force->inumeric(FLERR,strtok(NULL," \t\n\r\f"));
for (j = 0; j < ndouble; j++)
dvalues[j] = force->numeric(FLERR,strtok(NULL," \t\n\r\f"));
avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);
} else {
nvalues = ninteger + ndouble; // number of values to skip
for (j = 0; j < nvalues; j++)
strtok(NULL," \t\n\r\f");
}
}
delete [] ivalues;
delete [] dvalues;
}
/* ----------------------------------------------------------------------
allocate arrays of length ntypes
only done after ntypes is set
------------------------------------------------------------------------- */
void Atom::allocate_type_arrays()
{
if (avec->mass_type) {
mass = new double[ntypes+1];
mass_setflag = new int[ntypes+1];
for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
}
}
/* ----------------------------------------------------------------------
set a mass and flag it as set
called from reading of data file
type_offset may be used when reading multiple data files
------------------------------------------------------------------------- */
void Atom::set_mass(const char *str, int type_offset)
{
if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
int itype;
double mass_one;
int n = sscanf(str,"%d %lg",&itype,&mass_one);
if (n != 2) error->all(FLERR,"Invalid mass line in data file");
itype += type_offset;
if (itype < 1 || itype > ntypes)
error->all(FLERR,"Invalid type for mass set");
mass[itype] = mass_one;
mass_setflag[itype] = 1;
if (mass[itype] <= 0.0) error->all(FLERR,"Invalid mass value");
}
/* ----------------------------------------------------------------------
set a mass and flag it as set
called from EAM pair routine
------------------------------------------------------------------------- */
void Atom::set_mass(int itype, double value)
{
if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
if (itype < 1 || itype > ntypes)
error->all(FLERR,"Invalid type for mass set");
mass[itype] = value;
mass_setflag[itype] = 1;
if (mass[itype] <= 0.0) error->all(FLERR,"Invalid mass value");
}
/* ----------------------------------------------------------------------
set one or more masses and flag them as set
called from reading of input script
------------------------------------------------------------------------- */
void Atom::set_mass(int narg, char **arg)
{
if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
int lo,hi;
force->bounds(arg[0],ntypes,lo,hi);
if (lo < 1 || hi > ntypes) error->all(FLERR,"Invalid type for mass set");
for (int itype = lo; itype <= hi; itype++) {
mass[itype] = atof(arg[1]);
mass_setflag[itype] = 1;
if (mass[itype] <= 0.0) error->all(FLERR,"Invalid mass value");
}
}
/* ----------------------------------------------------------------------
set all masses as read in from restart file
------------------------------------------------------------------------- */
void Atom::set_mass(double *values)
{
for (int itype = 1; itype <= ntypes; itype++) {
mass[itype] = values[itype];
mass_setflag[itype] = 1;
}
}
/* ----------------------------------------------------------------------
check that all masses have been set
------------------------------------------------------------------------- */
void Atom::check_mass()
{
if (mass == NULL) return;
for (int itype = 1; itype <= ntypes; itype++)
if (mass_setflag[itype] == 0) error->all(FLERR,"All masses are not set");
}
/* ----------------------------------------------------------------------
check that radii of all particles of itype are the same
return 1 if true, else return 0
also return the radius value for that type
------------------------------------------------------------------------- */
int Atom::radius_consistency(int itype, double &rad)
{
double value = -1.0;
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (type[i] != itype) continue;
if (value < 0.0) value = radius[i];
else if (value != radius[i]) flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall) return 0;
MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world);
return 1;
}
/* ----------------------------------------------------------------------
check that shape of all particles of itype are the same
return 1 if true, else return 0
also return the 3 shape params for itype
------------------------------------------------------------------------- */
int Atom::shape_consistency(int itype,
double &shapex, double &shapey, double &shapez)
{
double zero[3] = {0.0, 0.0, 0.0};
double one[3] = {-1.0, -1.0, -1.0};
double *shape;
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int flag = 0;
for (int i = 0; i < nlocal; i++) {
if (type[i] != itype) continue;
if (ellipsoid[i] < 0) shape = zero;
else shape = bonus[ellipsoid[i]].shape;
if (one[0] < 0.0) {
one[0] = shape[0];
one[1] = shape[1];
one[2] = shape[2];
} else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2])
flag = 1;
}
int flagall;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
if (flagall) return 0;
double oneall[3];
MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world);
shapex = oneall[0];
shapey = oneall[1];
shapez = oneall[2];
return 1;
}
/* ----------------------------------------------------------------------
add a new molecule template = set of molecules
------------------------------------------------------------------------- */
void Atom::add_molecule(int narg, char **arg)
{
if (narg < 1) error->all(FLERR,"Illegal molecule command");
if (find_molecule(arg[0]) >= 0)
error->all(FLERR,"Reuse of molecule template ID");
// 1st molecule in set stores nset = # of mols, others store nset = 0
// ifile = count of molecules in set
// index = argument index where next molecule starts, updated by constructor
int ifile = 1;
int index = 1;
while (1) {
molecules = (Molecule **)
memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *),
"atom::molecules");
molecules[nmolecule] = new Molecule(lmp,narg,arg,index);
molecules[nmolecule]->nset = 0;
molecules[nmolecule-ifile+1]->nset++;
nmolecule++;
if (molecules[nmolecule-1]->last) break;
ifile++;
}
}
/* ----------------------------------------------------------------------
find first molecule in set with template ID
return -1 if does not exist
------------------------------------------------------------------------- */
int Atom::find_molecule(char *id)
{
int imol;
for (imol = 0; imol < nmolecule; imol++)
if (strcmp(id,molecules[imol]->id) == 0) return imol;
return -1;
}
/* ----------------------------------------------------------------------
add info to current atom ilocal from molecule template onemol and its iatom
offset = atom ID preceeding IDs of atoms in this molecule
called by fixes and commands that add molecules
------------------------------------------------------------------------- */
void Atom::add_molecule_atom(Molecule *onemol, int iatom,
int ilocal, tagint offset)
{
if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
else if (rmass_flag)
rmass[ilocal] = 4.0*MY_PI/3.0 *
radius[ilocal]*radius[ilocal]*radius[ilocal];
if (onemol->bodyflag) {
body[ilocal] = 0; // as if a body read from data file
onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody,
onemol->ibodyparams,onemol->dbodyparams);
onemol->avec_body->set_quat(ilocal,onemol->quat_external);
}
if (molecular != 1) return;
// add bond topology info
// for molecular atom styles, but not atom style template
if (avec->bonds_allow) {
num_bond[ilocal] = onemol->num_bond[iatom];
for (int i = 0; i < num_bond[ilocal]; i++) {
bond_type[ilocal][i] = onemol->bond_type[iatom][i];
bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset;
}
}
if (avec->angles_allow) {
num_angle[ilocal] = onemol->num_angle[iatom];
for (int i = 0; i < num_angle[ilocal]; i++) {
angle_type[ilocal][i] = onemol->angle_type[iatom][i];
angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset;
angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset;
angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset;
}
}
if (avec->dihedrals_allow) {
num_dihedral[ilocal] = onemol->num_dihedral[iatom];
for (int i = 0; i < num_dihedral[ilocal]; i++) {
dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i];
dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset;
dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset;
dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset;
dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset;
}
}
if (avec->impropers_allow) {
num_improper[ilocal] = onemol->num_improper[iatom];
for (int i = 0; i < num_improper[ilocal]; i++) {
improper_type[ilocal][i] = onemol->improper_type[iatom][i];
improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset;
improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset;
improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset;
improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset;
}
}
if (onemol->specialflag) {
nspecial[ilocal][0] = onemol->nspecial[iatom][0];
nspecial[ilocal][1] = onemol->nspecial[iatom][1];
int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2];
for (int i = 0; i < n; i++)
special[ilocal][i] = onemol->special[iatom][i] + offset;
}
}
/* ----------------------------------------------------------------------
reorder owned atoms so those in firstgroup appear first
called by comm->exchange() if atom_modify first group is set
only owned atoms exist at this point, no ghost atoms
------------------------------------------------------------------------- */
void Atom::first_reorder()
{
// insure there is one extra atom location at end of arrays for swaps
if (nlocal == nmax) avec->grow(0);
// loop over owned atoms
// nfirst = index of first atom not in firstgroup
// when find firstgroup atom out of place, swap it with atom nfirst
int bitmask = group->bitmask[firstgroup];
nfirst = 0;
while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & bitmask && i > nfirst) {
avec->copy(i,nlocal,0);
avec->copy(nfirst,i,0);
avec->copy(nlocal,nfirst,0);
while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
}
}
}
/* ----------------------------------------------------------------------
perform spatial sort of atoms within my sub-domain
always called between comm->exchange() and comm->borders()
don't have to worry about clearing/setting atom->map since done in comm
------------------------------------------------------------------------- */
void Atom::sort()
{
int i,m,n,ix,iy,iz,ibin,empty;
// set next timestep for sorting to take place
nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;
// download data from GPU if necessary
if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->downloadAll();
// re-setup sort bins if needed
if (domain->box_change) setup_sort_bins();
if (nbins == 1) return;
// reallocate per-atom vectors if needed
if (nlocal > maxnext) {
memory->destroy(next);
memory->destroy(permute);
maxnext = atom->nmax;
memory->create(next,maxnext,"atom:next");
memory->create(permute,maxnext,"atom:permute");
}
// insure there is one extra atom location at end of arrays for swaps
if (nlocal == nmax) avec->grow(0);
// bin atoms in reverse order so linked list will be in forward order
for (i = 0; i < nbins; i++) binhead[i] = -1;
for (i = nlocal-1; i >= 0; i--) {
ix = static_cast<int> ((x[i][0]-bboxlo[0])*bininvx);
iy = static_cast<int> ((x[i][1]-bboxlo[1])*bininvy);
iz = static_cast<int> ((x[i][2]-bboxlo[2])*bininvz);
ix = MAX(ix,0);
iy = MAX(iy,0);
iz = MAX(iz,0);
ix = MIN(ix,nbinx-1);
iy = MIN(iy,nbiny-1);
iz = MIN(iz,nbinz-1);
ibin = iz*nbiny*nbinx + iy*nbinx + ix;
next[i] = binhead[ibin];
binhead[ibin] = i;
}
// permute = desired permutation of atoms
// permute[I] = J means Ith new atom will be Jth old atom
n = 0;
for (m = 0; m < nbins; m++) {
i = binhead[m];
while (i >= 0) {
permute[n++] = i;
i = next[i];
}
}
// current = current permutation, just reuse next vector
// current[I] = J means Ith current atom is Jth old atom
int *current = next;
for (i = 0; i < nlocal; i++) current[i] = i;
// reorder local atom list, when done, current = permute
// perform "in place" using copy() to extra atom location at end of list
// inner while loop processes one cycle of the permutation
// copy before inner-loop moves an atom to end of atom list
// copy after inner-loop moves atom at end of list back into list
// empty = location in atom list that is currently empty
for (i = 0; i < nlocal; i++) {
if (current[i] == permute[i]) continue;
avec->copy(i,nlocal,0);
empty = i;
while (permute[empty] != i) {
avec->copy(permute[empty],empty,0);
empty = current[empty] = permute[empty];
}
avec->copy(nlocal,empty,0);
current[empty] = permute[empty];
}
// upload data back to GPU if necessary
if (lmp->cuda && !lmp->cuda->oncpu) lmp->cuda->uploadAll();
// sanity check that current = permute
//int flag = 0;
//for (i = 0; i < nlocal; i++)
// if (current[i] != permute[i]) flag = 1;
//int flagall;
//MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
//if (flagall) error->all(FLERR,"Atom sort did not operate correctly");
}
/* ----------------------------------------------------------------------
setup bins for spatial sorting of atoms
------------------------------------------------------------------------- */
void Atom::setup_sort_bins()
{
// binsize:
// user setting if explicitly set
// 1/2 of neighbor cutoff for non-CUDA
// CUDA_CHUNK atoms/proc for CUDA
// check if neighbor cutoff = 0.0
double binsize;
if (userbinsize > 0.0) binsize = userbinsize;
else if (!lmp->cuda) binsize = 0.5 * neighbor->cutneighmax;
else {
if (domain->dimension == 3) {
double vol = (domain->boxhi[0]-domain->boxlo[0]) *
(domain->boxhi[1]-domain->boxlo[1]) *
(domain->boxhi[2]-domain->boxlo[2]);
binsize = pow(1.0*CUDA_CHUNK/natoms*vol,1.0/3.0);
} else {
double area = (domain->boxhi[0]-domain->boxlo[0]) *
(domain->boxhi[1]-domain->boxlo[1]);
binsize = pow(1.0*CUDA_CHUNK/natoms*area,1.0/2.0);
}
}
if (binsize == 0.0) error->all(FLERR,"Atom sorting has bin size = 0.0");
double bininv = 1.0/binsize;
// nbin xyz = local bins
// bbox lo/hi = bounding box of my sub-domain
if (domain->triclinic)
domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
else {
bboxlo[0] = domain->sublo[0];
bboxlo[1] = domain->sublo[1];
bboxlo[2] = domain->sublo[2];
bboxhi[0] = domain->subhi[0];
bboxhi[1] = domain->subhi[1];
bboxhi[2] = domain->subhi[2];
}
nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) * bininv);
nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) * bininv);
nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) * bininv);
if (domain->dimension == 2) nbinz = 1;
if (nbinx == 0) nbinx = 1;
if (nbiny == 0) nbiny = 1;
if (nbinz == 0) nbinz = 1;
bininvx = nbinx / (bboxhi[0]-bboxlo[0]);
bininvy = nbiny / (bboxhi[1]-bboxlo[1]);
bininvz = nbinz / (bboxhi[2]-bboxlo[2]);
if (1.0*nbinx*nbiny*nbinz > INT_MAX)
error->one(FLERR,"Too many atom sorting bins");
nbins = nbinx*nbiny*nbinz;
// reallocate per-bin memory if needed
if (nbins > maxbin) {
memory->destroy(binhead);
maxbin = nbins;
memory->create(binhead,maxbin,"atom:binhead");
}
}
/* ----------------------------------------------------------------------
register a callback to a fix so it can manage atom-based arrays
happens when fix is created
flag = 0 for grow, 1 for restart, 2 for border comm
------------------------------------------------------------------------- */
void Atom::add_callback(int flag)
{
int ifix;
// find the fix
// if find NULL ptr:
// it's this one, since it is being replaced and has just been deleted
// at this point in re-creation
// if don't find NULL ptr:
// i is set to nfix = new one currently being added at end of list
for (ifix = 0; ifix < modify->nfix; ifix++)
if (modify->fix[ifix] == NULL) break;
// add callback to lists, reallocating if necessary
if (flag == 0) {
if (nextra_grow == nextra_grow_max) {
nextra_grow_max += DELTA;
memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow");
}
extra_grow[nextra_grow] = ifix;
nextra_grow++;
} else if (flag == 1) {
if (nextra_restart == nextra_restart_max) {
nextra_restart_max += DELTA;
memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart");
}
extra_restart[nextra_restart] = ifix;
nextra_restart++;
} else if (flag == 2) {
if (nextra_border == nextra_border_max) {
nextra_border_max += DELTA;
memory->grow(extra_border,nextra_border_max,"atom:extra_border");
}
extra_border[nextra_border] = ifix;
nextra_border++;
}
}
/* ----------------------------------------------------------------------
unregister a callback to a fix
happens when fix is deleted, called by its destructor
flag = 0 for grow, 1 for restart
------------------------------------------------------------------------- */
void Atom::delete_callback(const char *id, int flag)
{
int ifix;
for (ifix = 0; ifix < modify->nfix; ifix++)
if (strcmp(id,modify->fix[ifix]->id) == 0) break;
// compact the list of callbacks
if (flag == 0) {
int match;
for (match = 0; match < nextra_grow; match++)
if (extra_grow[match] == ifix) break;
for (int i = match; i < nextra_grow-1; i++)
extra_grow[i] = extra_grow[i+1];
nextra_grow--;
} else if (flag == 1) {
int match;
for (match = 0; match < nextra_restart; match++)
if (extra_restart[match] == ifix) break;
for (int i = match; i < nextra_restart-1; i++)
extra_restart[i] = extra_restart[i+1];
nextra_restart--;
} else if (flag == 2) {
int match;
for (match = 0; match < nextra_border; match++)
if (extra_border[match] == ifix) break;
for (int i = match; i < nextra_border-1; i++)
extra_border[i] = extra_border[i+1];
nextra_border--;
}
}
/* ----------------------------------------------------------------------
decrement ptrs in callback lists to fixes beyond the deleted ifix
happens after fix is deleted
------------------------------------------------------------------------- */
void Atom::update_callback(int ifix)
{
for (int i = 0; i < nextra_grow; i++)
if (extra_grow[i] > ifix) extra_grow[i]--;
for (int i = 0; i < nextra_restart; i++)
if (extra_restart[i] > ifix) extra_restart[i]--;
for (int i = 0; i < nextra_border; i++)
if (extra_border[i] > ifix) extra_border[i]--;
}
/* ----------------------------------------------------------------------
find custom per-atom vector with name
return index if found, and flag = 0/1 for int/double
return -1 if not found
------------------------------------------------------------------------- */
int Atom::find_custom(char *name, int &flag)
{
for (int i = 0; i < nivector; i++)
if (iname[i] && strcmp(iname[i],name) == 0) {
flag = 0;
return i;
}
for (int i = 0; i < ndvector; i++)
if (dname[i] && strcmp(dname[i],name) == 0) {
flag = 1;
return i;
}
return -1;
}
/* ----------------------------------------------------------------------
add a custom variable with name of type flag = 0/1 for int/double
assumes name does not already exist
return index in ivector or dvector of its location
------------------------------------------------------------------------- */
int Atom::add_custom(char *name, int flag)
{
int index;
if (flag == 0) {
index = nivector;
nivector++;
iname = (char **) memory->srealloc(iname,nivector*sizeof(char *),
"atom:iname");
int n = strlen(name) + 1;
iname[index] = new char[n];
strcpy(iname[index],name);
ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),
"atom:ivector");
memory->create(ivector[index],nmax,"atom:ivector");
} else {
index = ndvector;
ndvector++;
dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *),
"atom:dname");
int n = strlen(name) + 1;
dname[index] = new char[n];
strcpy(dname[index],name);
dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),
"atom:dvector");
memory->create(dvector[index],nmax,"atom:dvector");
}
return index;
}
/* ----------------------------------------------------------------------
remove a custom variable of type flag = 0/1 for int/double at index
free memory for vector and name and set ptrs to NULL
ivector/dvector and iname/dname lists never shrink
------------------------------------------------------------------------- */
void Atom::remove_custom(int flag, int index)
{
if (flag == 0) {
memory->destroy(ivector[index]);
ivector[index] = NULL;
delete [] iname[index];
iname[index] = NULL;
} else {
memory->destroy(dvector[index]);
dvector[index] = NULL;
delete [] dname[index];
dname[index] = NULL;
}
}
/* ----------------------------------------------------------------------
return a pointer to a named internal variable
if don't recognize name, return NULL
customize by adding names
------------------------------------------------------------------------- */
void *Atom::extract(char *name)
{
if (strcmp(name,"mass") == 0) return (void *) mass;
if (strcmp(name,"id") == 0) return (void *) tag;
if (strcmp(name,"type") == 0) return (void *) type;
if (strcmp(name,"mask") == 0) return (void *) mask;
if (strcmp(name,"image") == 0) return (void *) image;
if (strcmp(name,"x") == 0) return (void *) x;
if (strcmp(name,"v") == 0) return (void *) v;
if (strcmp(name,"f") == 0) return (void *) f;
if (strcmp(name,"molecule") == 0) return (void *) molecule;
if (strcmp(name,"q") == 0) return (void *) q;
if (strcmp(name,"mu") == 0) return (void *) mu;
if (strcmp(name,"omega") == 0) return (void *) omega;
if (strcmp(name,"angmom") == 0) return (void *) angmom;
if (strcmp(name,"torque") == 0) return (void *) torque;
if (strcmp(name,"radius") == 0) return (void *) radius;
if (strcmp(name,"rmass") == 0) return (void *) rmass;
if (strcmp(name,"ellipsoid") == 0) return (void *) ellipsoid;
if (strcmp(name,"line") == 0) return (void *) line;
if (strcmp(name,"tri") == 0) return (void *) tri;
if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
if (strcmp(name,"s0") == 0) return (void *) s0;
if (strcmp(name,"x0") == 0) return (void *) x0;
if (strcmp(name,"spin") == 0) return (void *) spin;
if (strcmp(name,"eradius") == 0) return (void *) eradius;
if (strcmp(name,"ervel") == 0) return (void *) ervel;
if (strcmp(name,"erforce") == 0) return (void *) erforce;
if (strcmp(name,"ervelforce") == 0) return (void *) ervelforce;
if (strcmp(name,"cs") == 0) return (void *) cs;
if (strcmp(name,"csforce") == 0) return (void *) csforce;
if (strcmp(name,"vforce") == 0) return (void *) vforce;
if (strcmp(name,"etag") == 0) return (void *) etag;
if (strcmp(name,"rho") == 0) return (void *) rho;
if (strcmp(name,"drho") == 0) return (void *) drho;
if (strcmp(name,"e") == 0) return (void *) e;
if (strcmp(name,"de") == 0) return (void *) de;
if (strcmp(name,"cv") == 0) return (void *) cv;
if (strcmp(name,"vest") == 0) return (void *) vest;
if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius;
if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9;
if (strcmp(name, "smd_stress") == 0) return (void *) smd_stress;
if (strcmp(name, "eff_plastic_strain") == 0)
return (void *) eff_plastic_strain;
if (strcmp(name, "eff_plastic_strain_rate") == 0)
return (void *) eff_plastic_strain_rate;
if (strcmp(name, "damage") == 0) return (void *) damage;
+ if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
+
return NULL;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
call to avec tallies per-atom vectors
add in global to local mapping storage
------------------------------------------------------------------------- */
bigint Atom::memory_usage()
{
memlength = DELTA_MEMSTR;
memory->create(memstr,memlength,"atom:memstr");
memstr[0] = '\0';
bigint bytes = avec->memory_usage();
memory->destroy(memstr);
bytes += max_same*sizeof(int);
if (map_style == 1)
bytes += memory->usage(map_array,map_maxarray);
else if (map_style == 2) {
bytes += map_nbucket*sizeof(int);
bytes += map_nhash*sizeof(HashElem);
}
if (maxnext) {
bytes += memory->usage(next,maxnext);
bytes += memory->usage(permute,maxnext);
}
return bytes;
}
/* ----------------------------------------------------------------------
accumulate per-atom vec names in memstr, padded by spaces
return 1 if padded str is not already in memlist, else 0
------------------------------------------------------------------------- */
int Atom::memcheck(const char *str)
{
int n = strlen(str) + 3;
char *padded = new char[n];
strcpy(padded," ");
strcat(padded,str);
strcat(padded," ");
if (strstr(memstr,padded)) {
delete [] padded;
return 0;
}
if (strlen(memstr) + n >= memlength) {
memlength += DELTA_MEMSTR;
memory->grow(memstr,memlength,"atom:memstr");
}
strcat(memstr,padded);
delete [] padded;
return 1;
}
diff --git a/src/atom.h b/src/atom.h
index c0b69c952..39611adf9 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -1,493 +1,500 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifndef LMP_ATOM_H
#define LMP_ATOM_H
#include "pointers.h"
namespace LAMMPS_NS {
class Atom : protected Pointers {
public:
char *atom_style;
class AtomVec *avec;
// atom counts
bigint natoms; // total # of atoms in system, could be 0
// natoms may not be current if atoms lost
int nlocal,nghost; // # of owned and ghost atoms on this proc
int nmax; // max # of owned+ghost in arrays on this proc
int tag_enable; // 0/1 if atom ID tags are defined
int molecular; // 0 = atomic, 1 = standard molecular system,
// 2 = molecule template system
bigint nbonds,nangles,ndihedrals,nimpropers;
int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
int extra_bond_per_atom,extra_angle_per_atom;
int extra_dihedral_per_atom,extra_improper_per_atom;
int firstgroup; // store atoms in this group first, -1 if unset
int nfirst; // # of atoms in first group on this proc
char *firstgroupname; // group-ID to store first, NULL if unset
// per-atom arrays
// customize by adding new array
tagint *tag;
int *type,*mask;
imageint *image;
double **x,**v,**f;
tagint *molecule;
int *molindex,*molatom;
double *q,**mu;
double **omega,**angmom,**torque;
double *radius,*rmass;
int *ellipsoid,*line,*tri,*body;
// PERI package
double *vfrac,*s0;
double **x0;
// USER-EFF and USER-AWPMD packages
int *spin;
double *eradius,*ervel,*erforce,*ervelforce;
double *cs,*csforce,*vforce;
int *etag;
// USER-SPH package
double *rho,*drho,*e,*de,*cv;
double **vest;
// USER-SMD package
double *contact_radius;
double **smd_data_9;
double **smd_stress;
double *eff_plastic_strain;
double *eff_plastic_strain_rate;
double *damage;
+ // USER-DPD package
+
+ double *uCond, *uMech, *uChem, *uCGnew, *uCG;
+ double *duCond, *duMech, *duChem;
+ double *dpdTheta;
+
// molecular info
int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
tagint **special; // IDs of 1-2,1-3,1-4 neighs of each atom
int maxspecial; // special[nlocal][maxspecial]
int *num_bond;
int **bond_type;
tagint **bond_atom;
int *num_angle;
int **angle_type;
tagint **angle_atom1,**angle_atom2,**angle_atom3;
int *num_dihedral;
int **dihedral_type;
tagint **dihedral_atom1,**dihedral_atom2,**dihedral_atom3,**dihedral_atom4;
int *num_improper;
int **improper_type;
tagint **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
// custom arrays used by fix property/atom
int **ivector;
double **dvector;
char **iname,**dname;
int nivector,ndvector;
// used by USER-CUDA to flag used per-atom arrays
unsigned int datamask;
unsigned int datamask_ext;
// atom style and per-atom array existence flags
// customize by adding new flag
int sphere_flag,ellipsoid_flag,line_flag,tri_flag,body_flag;
int peri_flag,electron_flag;
int ecp_flag;
int wavepacket_flag,sph_flag;
int molecule_flag,molindex_flag,molatom_flag;
int q_flag,mu_flag;
int rmass_flag,radius_flag,omega_flag,torque_flag,angmom_flag;
int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag;
int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag;
int rho_flag,e_flag,cv_flag,vest_flag;
+ int dpd_flag;
// USER-SMD package
int smd_flag;
int contact_radius_flag;
int smd_data_9_flag;
int smd_stress_flag;
int x0_flag;
int eff_plastic_strain_flag;
int eff_plastic_strain_rate_flag;
int damage_flag;
// Peridynamics scale factor, used by dump cfg
double pdscale;
// molecule templates
// each template can be a set of consecutive molecules
// each with same ID (stored in molecules)
// 1st molecule in template stores nset = # in set
int nmolecule;
class Molecule **molecules;
// extra peratom info in restart file destined for fix & diag
double **extra;
// per-type arrays
double *mass;
int *mass_setflag;
// callback ptrs for atom arrays managed by fix classes
int nextra_grow,nextra_restart,nextra_border; // # of callbacks of each type
int *extra_grow,*extra_restart,*extra_border; // index of fix to callback to
int nextra_grow_max,nextra_restart_max; // size of callback lists
int nextra_border_max;
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
tagint map_tag_max; // max atom ID that map() is setup for
// spatial sorting of atoms
int sortfreq; // sort atoms every this many steps, 0 = off
bigint nextsort; // next timestep to sort on
double userbinsize; // requested sort bin size
// indices of atoms with same ID
int *sametag; // sametag[I] = next atom with same ID, -1 if no more
// functions
Atom(class LAMMPS *);
~Atom();
void settings(class Atom *);
void create_avec(const char *, int, char **, int);
virtual class AtomVec *new_avec(const char *, int, int &);
void init();
void setup();
class AtomVec *style_match(const char *);
void modify_params(int, char **);
void tag_check();
void tag_extend();
int tag_consecutive();
int parse_data(const char *);
int count_words(const char *);
int count_words(const char *, char *);
void deallocate_topology();
void data_atoms(int, char *, tagint, int, int, double *);
void data_vels(int, char *, tagint);
void data_bonds(int, char *, int *, tagint, int);
void data_angles(int, char *, int *, tagint, int);
void data_dihedrals(int, char *, int *, tagint, int);
void data_impropers(int, char *, int *, tagint, int);
void data_bonus(int, char *, class AtomVec *, tagint);
void data_bodies(int, char *, class AtomVecBody *, tagint);
virtual void allocate_type_arrays();
void set_mass(const char *, int);
void set_mass(int, double);
void set_mass(int, char **);
void set_mass(double *);
void check_mass();
int radius_consistency(int, double &);
int shape_consistency(int, double &, double &, double &);
void add_molecule(int, char **);
int find_molecule(char *);
void add_molecule_atom(class Molecule *, int, int, tagint);
void first_reorder();
virtual void sort();
void add_callback(int);
void delete_callback(const char *, int);
void update_callback(int);
int find_custom(char *, int &);
int add_custom(char *, int);
void remove_custom(int, int);
virtual void sync_modify(ExecutionSpace, unsigned int, unsigned int) {}
void *extract(char *);
inline int* get_map_array() {return map_array;};
inline int get_map_size() {return map_tag_max+1;};
bigint memory_usage();
int memcheck(const char *);
// functions for global to local ID mapping
// map lookup function inlined for efficiency
// return -1 if no map defined
inline int map(tagint global) {
if (map_style == 1) return map_array[global];
else if (map_style == 2) return map_find_hash(global);
else return -1;
};
void map_init(int check = 1);
void map_clear();
void map_set();
void map_one(tagint, int);
int map_style_set();
void map_delete();
int map_find_hash(tagint);
protected:
// global to local ID mapping
int *map_array; // direct map via array that holds map_tag_max
int map_maxarray; // allocated size of map_array (1 larger than this)
struct HashElem { // hashed map
tagint global; // key to search on = global ID
int local; // value associated with key = local index
int next; // next entry in this bucket, -1 if last
};
int map_nhash; // # of entries hash table can hold
int map_nused; // # of actual entries in hash table
int map_free; // ptr to 1st unused entry in hash table
int map_nbucket; // # of hash buckets
int *map_bucket; // ptr to 1st entry in each bucket
HashElem *map_hash; // hash table
int max_same; // allocated size of sametag
// spatial sorting of atoms
int nbins; // # of sorting bins
int nbinx,nbiny,nbinz; // bins in each dimension
int maxbin; // max # of bins
int maxnext; // max size of next,permute
int *binhead; // 1st atom in each bin
int *next; // next atom in bin
int *permute; // permutation vector
double bininvx,bininvy,bininvz; // inverse actual bin sizes
double bboxlo[3],bboxhi[3]; // bounding box of my sub-domain
int memlength; // allocated size of memstr
char *memstr; // string of array names already counted
void setup_sort_bins();
int next_prime(int);
};
}
#endif
/* ERROR/WARNING messages:
E: Atom IDs must be used for molecular systems
Atom IDs are used to identify and find partner atoms in bonds.
E: Unknown atom style
The choice of atom style is unknown.
E: Could not find atom_modify first group ID
Self-explanatory.
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Atom_modify id command after simulation box is defined
The atom_modify id command cannot be used after a read_data,
read_restart, or create_box command.
E: Atom_modify map command after simulation box is defined
The atom_modify map command cannot be used after a read_data,
read_restart, or create_box command.
E: Atom_modify sort and first options cannot be used together
Self-explanatory.
E: Atom ID is negative
Self-explanatory.
E: Atom ID is too big
The limit on atom IDs is set by the SMALLBIG, BIGBIG, SMALLSMALL
setting in your Makefile. See Section_start 2.2 of the manual for
more details.
E: Atom ID is zero
Either all atoms IDs must be zero or none of them.
E: Not all atom IDs are 0
Either all atoms IDs must be zero or none of them.
E: New atom IDs exceed maximum allowed ID
See the setting for tagint in the src/lmptype.h file.
E: Incorrect atom format in data file
Number of values per atom line in the data file is not consistent with
the atom style.
E: Incorrect velocity format in data file
Each atom style defines a format for the Velocity section
of the data file. The read-in lines do not match.
E: Invalid atom ID in Velocities section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Invalid atom ID in Bonds section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Invalid bond type in Bonds section of data file
Bond type must be positive integer and within range of specified bond
types.
E: Invalid atom ID in Angles section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Invalid angle type in Angles section of data file
Angle type must be positive integer and within range of specified angle
types.
E: Invalid atom ID in Dihedrals section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Invalid dihedral type in Dihedrals section of data file
Dihedral type must be positive integer and within range of specified
dihedral types.
E: Invalid atom ID in Impropers section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Invalid improper type in Impropers section of data file
Improper type must be positive integer and within range of specified
improper types.
E: Incorrect bonus data format in data file
See the read_data doc page for a description of how various kinds of
bonus data must be formatted for certain atom styles.
E: Invalid atom ID in Bonus section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Invalid atom ID in Bodies section of data file
Atom IDs must be positive integers and within range of defined
atoms.
E: Cannot set mass for this atom style
This atom style does not support mass settings for each atom type.
Instead they are defined on a per-atom basis in the data file.
E: Invalid mass line in data file
Self-explanatory.
E: Invalid type for mass set
Mass command must set a type from 1-N where N is the number of atom
types.
E: Invalid mass value
Self-explanatory.
E: All masses are not set
For atom styles that define masses for each atom type, all masses must
be set in the data file or by the mass command before running a
simulation. They must also be set before using the velocity
command.
E: Reuse of molecule template ID
The template IDs must be unique.
E: Atom sort did not operate correctly
This is an internal LAMMPS error. Please report it to the
developers.
E: Atom sorting has bin size = 0.0
The neighbor cutoff is being used as the bin size, but it is zero.
Thus you must explicitly list a bin size in the atom_modify sort
command or turn off sorting.
E: Too many atom sorting bins
This is likely due to an immense simulation box that has blown up
to a large size.
*/
diff --git a/src/library.cpp b/src/library.cpp
index 5801cedc5..44734f95f 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -1,557 +1,559 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
// C or Fortran style library interface to LAMMPS
// customize by adding new LAMMPS-specific functions
#include <mpi.h>
#include <string.h>
#include <stdlib.h>
#include "library.h"
#include "lammps.h"
#include "universe.h"
#include "input.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
#include "group.h"
#include "input.h"
#include "variable.h"
#include "modify.h"
#include "compute.h"
#include "fix.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ----------------------------------------------------------------------
create an instance of LAMMPS and return pointer to it
pass in command-line args and MPI communicator to run on
------------------------------------------------------------------------- */
void lammps_open(int argc, char **argv, MPI_Comm communicator, void **ptr)
{
LAMMPS *lmp = new LAMMPS(argc,argv,communicator);
*ptr = (void *) lmp;
}
/* ----------------------------------------------------------------------
create an instance of LAMMPS and return pointer to it
caller doesn't know MPI communicator, so use MPI_COMM_WORLD
intialize MPI if needed
------------------------------------------------------------------------- */
void lammps_open_no_mpi(int argc, char **argv, void **ptr)
{
int flag;
MPI_Initialized(&flag);
if (!flag) {
int argc = 0;
char **argv = NULL;
MPI_Init(&argc,&argv);
}
MPI_Comm communicator = MPI_COMM_WORLD;
LAMMPS *lmp = new LAMMPS(argc,argv,communicator);
*ptr = (void *) lmp;
}
/* ----------------------------------------------------------------------
destruct an instance of LAMMPS
------------------------------------------------------------------------- */
void lammps_close(void *ptr)
{
LAMMPS *lmp = (LAMMPS *) ptr;
delete lmp;
}
/* ----------------------------------------------------------------------
get the numerical representation of the current LAMMPS version
------------------------------------------------------------------------- */
int lammps_version(void *ptr)
{
LAMMPS *lmp = (LAMMPS *) ptr;
return atoi(lmp->universe->num_ver);
}
/* ----------------------------------------------------------------------
process an input script in filename str
------------------------------------------------------------------------- */
void lammps_file(void *ptr, char *str)
{
LAMMPS *lmp = (LAMMPS *) ptr;
lmp->input->file(str);
}
/* ----------------------------------------------------------------------
process a single input command in str
------------------------------------------------------------------------- */
char *lammps_command(void *ptr, char *str)
{
LAMMPS *lmp = (LAMMPS *) ptr;
return lmp->input->one(str);
}
/* ----------------------------------------------------------------------
clean-up function to free memory allocated by lib and returned to caller
------------------------------------------------------------------------- */
void lammps_free(void *ptr)
{
free(ptr);
}
/* ----------------------------------------------------------------------
add LAMMPS-specific library functions
all must receive LAMMPS pointer as argument
customize by adding a function here and in library.h header file
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
extract a pointer to an internal LAMMPS global entity
name = desired quantity, e.g. dt or boxyhi or natoms
returns a void pointer to the entity
which the caller can cast to the proper data type
returns a NULL if name not listed below
customize by adding names
------------------------------------------------------------------------- */
void *lammps_extract_global(void *ptr, char *name)
{
LAMMPS *lmp = (LAMMPS *) ptr;
if (strcmp(name,"dt") == 0) return (void *) &lmp->update->dt;
if (strcmp(name,"boxxlo") == 0) return (void *) &lmp->domain->boxlo[0];
if (strcmp(name,"boxxhi") == 0) return (void *) &lmp->domain->boxhi[0];
if (strcmp(name,"boxylo") == 0) return (void *) &lmp->domain->boxlo[1];
if (strcmp(name,"boxyhi") == 0) return (void *) &lmp->domain->boxhi[1];
if (strcmp(name,"boxzlo") == 0) return (void *) &lmp->domain->boxlo[2];
if (strcmp(name,"boxzhi") == 0) return (void *) &lmp->domain->boxhi[2];
if (strcmp(name,"xy") == 0) return (void *) &lmp->domain->xy;
if (strcmp(name,"xz") == 0) return (void *) &lmp->domain->xz;
if (strcmp(name,"yz") == 0) return (void *) &lmp->domain->yz;
if (strcmp(name,"natoms") == 0) return (void *) &lmp->atom->natoms;
if (strcmp(name,"nlocal") == 0) return (void *) &lmp->atom->nlocal;
+ if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep;
+
return NULL;
}
/* ----------------------------------------------------------------------
extract a pointer to an internal LAMMPS atom-based entity
name = desired quantity, e.g. x or mass
returns a void pointer to the entity
which the caller can cast to the proper data type
returns a NULL if Atom::extract() does not recognize the name
customize by adding names to Atom::extract()
------------------------------------------------------------------------- */
void *lammps_extract_atom(void *ptr, char *name)
{
LAMMPS *lmp = (LAMMPS *) ptr;
return lmp->atom->extract(name);
}
/* ----------------------------------------------------------------------
extract a pointer to an internal LAMMPS compute-based entity
id = compute ID
style = 0 for global data, 1 for per-atom data, 2 for local data
type = 0 for scalar, 1 for vector, 2 for array
for global data, returns a pointer to the
compute's internal data structure for the entity
caller should cast it to (double *) for a scalar or vector
caller should cast it to (double **) for an array
for per-atom or local data, returns a pointer to the
compute's internal data structure for the entity
caller should cast it to (double *) for a vector
caller should cast it to (double **) for an array
returns a void pointer to the compute's internal data structure
for the entity which the caller can cast to the proper data type
returns a NULL if id is not recognized or style/type not supported
IMPORTANT: if the compute is not current it will be invoked
LAMMPS cannot easily check here if it is valid to invoke the compute,
so caller must insure that it is OK
------------------------------------------------------------------------- */
void *lammps_extract_compute(void *ptr, char *id, int style, int type)
{
LAMMPS *lmp = (LAMMPS *) ptr;
int icompute = lmp->modify->find_compute(id);
if (icompute < 0) return NULL;
Compute *compute = lmp->modify->compute[icompute];
if (style == 0) {
if (type == 0) {
if (!compute->scalar_flag) return NULL;
if (compute->invoked_scalar != lmp->update->ntimestep)
compute->compute_scalar();
return (void *) &compute->scalar;
}
if (type == 1) {
if (!compute->vector_flag) return NULL;
if (compute->invoked_vector != lmp->update->ntimestep)
compute->compute_vector();
return (void *) compute->vector;
}
if (type == 2) {
if (!compute->array_flag) return NULL;
if (compute->invoked_array != lmp->update->ntimestep)
compute->compute_array();
return (void *) compute->array;
}
}
if (style == 1) {
if (!compute->peratom_flag) return NULL;
if (type == 1) {
if (compute->invoked_peratom != lmp->update->ntimestep)
compute->compute_peratom();
return (void *) compute->vector_atom;
}
if (type == 2) {
if (compute->invoked_peratom != lmp->update->ntimestep)
compute->compute_peratom();
return (void *) compute->array_atom;
}
}
if (style == 2) {
if (!compute->local_flag) return NULL;
if (type == 1) {
if (compute->invoked_local != lmp->update->ntimestep)
compute->compute_local();
return (void *) compute->vector_local;
}
if (type == 2) {
if (compute->invoked_local != lmp->update->ntimestep)
compute->compute_local();
return (void *) compute->array_local;
}
}
return NULL;
}
/* ----------------------------------------------------------------------
extract a pointer to an internal LAMMPS fix-based entity
id = fix ID
style = 0 for global data, 1 for per-atom data, 2 for local data
type = 0 for scalar, 1 for vector, 2 for array
i,j = indices needed only to specify which global vector or array value
for global data, returns a pointer to a memory location
which is allocated by this function
which the caller can cast to a (double *) which points to the value
for per-atom or local data, returns a pointer to the
fix's internal data structure for the entity
caller should cast it to (double *) for a vector
caller should cast it to (double **) for an array
returns a NULL if id is not recognized or style/type not supported
IMPORTANT: for global data,
this function allocates a double to store the value in,
so the caller must free this memory to avoid a leak, e.g.
double *dptr = (double *) lammps_extract_fix();
double value = *dptr;
lammps_free(dptr);
IMPORTANT: LAMMPS cannot easily check here when info extracted from
the fix is valid, so caller must insure that it is OK
------------------------------------------------------------------------- */
void *lammps_extract_fix(void *ptr, char *id, int style, int type,
int i, int j)
{
LAMMPS *lmp = (LAMMPS *) ptr;
int ifix = lmp->modify->find_fix(id);
if (ifix < 0) return NULL;
Fix *fix = lmp->modify->fix[ifix];
if (style == 0) {
double *dptr = (double *) malloc(sizeof(double));
if (type == 0) {
if (!fix->scalar_flag) return NULL;
*dptr = fix->compute_scalar();
return (void *) dptr;
}
if (type == 1) {
if (!fix->vector_flag) return NULL;
*dptr = fix->compute_vector(i);
return (void *) dptr;
}
if (type == 2) {
if (!fix->array_flag) return NULL;
*dptr = fix->compute_array(i,j);
return (void *) dptr;
}
}
if (style == 1) {
if (!fix->peratom_flag) return NULL;
if (type == 1) return (void *) fix->vector_atom;
if (type == 2) return (void *) fix->array_atom;
}
if (style == 2) {
if (!fix->local_flag) return NULL;
if (type == 1) return (void *) fix->vector_local;
if (type == 2) return (void *) fix->array_local;
}
return NULL;
}
/* ----------------------------------------------------------------------
extract a pointer to an internal LAMMPS evaluated variable
name = variable name, must be equal-style or atom-style variable
group = group ID for evaluating an atom-style variable, else NULL
for equal-style variable, returns a pointer to a memory location
which is allocated by this function
which the caller can cast to a (double *) which points to the value
for atom-style variable, returns a pointer to the
vector of per-atom values on each processor,
which the caller can cast to a (double *) which points to the values
returns a NULL if name is not recognized or not equal-style or atom-style
IMPORTANT: for both equal-style and atom-style variables,
this function allocates memory to store the variable data in
so the caller must free this memory to avoid a leak
e.g. for equal-style variables
double *dptr = (double *) lammps_extract_variable();
double value = *dptr;
lammps_free(dptr);
e.g. for atom-style variables
double *vector = (double *) lammps_extract_variable();
use the vector values
lammps_free(vector);
IMPORTANT: LAMMPS cannot easily check here when it is valid to evaluate
the variable or any fixes or computes or thermodynamic info it references,
so caller must insure that it is OK
------------------------------------------------------------------------- */
void *lammps_extract_variable(void *ptr, char *name, char *group)
{
LAMMPS *lmp = (LAMMPS *) ptr;
int ivar = lmp->input->variable->find(name);
if (ivar < 0) return NULL;
if (lmp->input->variable->equalstyle(ivar)) {
double *dptr = (double *) malloc(sizeof(double));
*dptr = lmp->input->variable->compute_equal(ivar);
return (void *) dptr;
}
if (lmp->input->variable->atomstyle(ivar)) {
int igroup = lmp->group->find(group);
if (igroup < 0) return NULL;
int nlocal = lmp->atom->nlocal;
double *vector = (double *) malloc(nlocal*sizeof(double));
lmp->input->variable->compute_atom(ivar,igroup,vector,1,0);
return (void *) vector;
}
return NULL;
}
/* ----------------------------------------------------------------------
set the value of a STRING variable to str
return -1 if variable doesn't exist or not a STRING variable
return 0 for success
------------------------------------------------------------------------- */
int lammps_set_variable(void *ptr, char *name, char *str)
{
LAMMPS *lmp = (LAMMPS *) ptr;
int err = lmp->input->variable->set_string(name,str);
return err;
}
/* ----------------------------------------------------------------------
return the total number of atoms in the system
useful before call to lammps_get_atoms() so can pre-allocate vector
------------------------------------------------------------------------- */
int lammps_get_natoms(void *ptr)
{
LAMMPS *lmp = (LAMMPS *) ptr;
if (lmp->atom->natoms > MAXSMALLINT) return 0;
int natoms = static_cast<int> (lmp->atom->natoms);
return natoms;
}
/* ----------------------------------------------------------------------
gather the named atom-based entity across all processors
name = desired quantity, e.g. x or charge
type = 0 for integer values, 1 for double values
count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
return atom-based values in data, ordered by count, then by atom ID
e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
data must be pre-allocated by caller to correct length
------------------------------------------------------------------------- */
void lammps_gather_atoms(void *ptr, char *name,
int type, int count, void *data)
{
LAMMPS *lmp = (LAMMPS *) ptr;
// error if tags are not defined or not consecutive
int flag = 0;
if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
if (flag && lmp->comm->me == 0) {
lmp->error->warning(FLERR,"Library error in lammps_gather_atoms");
return;
}
int natoms = static_cast<int> (lmp->atom->natoms);
int i,j,offset;
void *vptr = lmp->atom->extract(name);
// copy = Natom length vector of per-atom values
// use atom ID to insert each atom's values into copy
// MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID
if (type == 0) {
int *vector = NULL;
int **array = NULL;
if (count == 1) vector = (int *) vptr;
else array = (int **) vptr;
int *copy = (int*) data;
for (i = 0; i < count*natoms; i++) copy[i] = 0;
tagint *tag = lmp->atom->tag;
int nlocal = lmp->atom->nlocal;
if (count == 1)
for (i = 0; i < nlocal; i++)
copy[tag[i]-1] = vector[i];
else
for (i = 0; i < nlocal; i++) {
offset = count*(tag[i]-1);
for (j = 0; j < count; j++)
copy[offset++] = array[i][0];
}
MPI_Allreduce(MPI_IN_PLACE,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
} else {
double *vector = NULL;
double **array = NULL;
if (count == 1) vector = (double *) vptr;
else array = (double **) vptr;
double *copy = (double*) data;
for (i = 0; i < count*natoms; i++) copy[i] = 0.0;
tagint *tag = lmp->atom->tag;
int nlocal = lmp->atom->nlocal;
if (count == 1) {
for (i = 0; i < nlocal; i++)
copy[tag[i]-1] = vector[i];
} else {
for (i = 0; i < nlocal; i++) {
offset = count*(tag[i]-1);
for (j = 0; j < count; j++)
copy[offset++] = array[i][j];
}
}
MPI_Allreduce(MPI_IN_PLACE,data,count*natoms,MPI_DOUBLE,MPI_SUM,lmp->world);
}
}
/* ----------------------------------------------------------------------
scatter the named atom-based entity across all processors
name = desired quantity, e.g. x or charge
type = 0 for integer values, 1 for double values
count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
data = atom-based values in data, ordered by count, then by atom ID
e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
------------------------------------------------------------------------- */
void lammps_scatter_atoms(void *ptr, char *name,
int type, int count, void *data)
{
LAMMPS *lmp = (LAMMPS *) ptr;
// error if tags are not defined or not consecutive or no atom map
int flag = 0;
if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
if (lmp->atom->map_style == 0) flag = 1;
if (flag && lmp->comm->me == 0) {
lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms");
return;
}
int natoms = static_cast<int> (lmp->atom->natoms);
int i,j,m,offset;
void *vptr = lmp->atom->extract(name);
// copy = Natom length vector of per-atom values
// use atom ID to insert each atom's values into copy
// MPI_Allreduce with MPI_SUM to merge into data, ordered by atom ID
if (type == 0) {
int *vector = NULL;
int **array = NULL;
if (count == 1) vector = (int *) vptr;
else array = (int **) vptr;
int *dptr = (int *) data;
if (count == 1) {
for (i = 0; i < natoms; i++)
if ((m = lmp->atom->map(i+1)) >= 0)
vector[m] = dptr[i];
} else {
for (i = 0; i < natoms; i++)
if ((m = lmp->atom->map(i+1)) >= 0) {
offset = count*i;
for (j = 0; j < count; j++)
array[m][j] = dptr[offset++];
}
}
} else {
double *vector = NULL;
double **array = NULL;
if (count == 1) vector = (double *) vptr;
else array = (double **) vptr;
double *dptr = (double *) data;
if (count == 1) {
for (i = 0; i < natoms; i++)
if ((m = lmp->atom->map(i+1)) >= 0)
vector[m] = dptr[i];
} else {
for (i = 0; i < natoms; i++) {
if ((m = lmp->atom->map(i+1)) >= 0) {
offset = count*i;
for (j = 0; j < count; j++)
array[m][j] = dptr[offset++];
}
}
}
}
}
diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp
index 7c9cd448c..9a7152a1b 100644
--- a/src/neigh_request.cpp
+++ b/src/neigh_request.cpp
@@ -1,217 +1,222 @@
/* ----------------------------------------------------------------------
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 "neigh_request.h"
#include "atom.h"
#include "memory.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
{
// default ID = 0
id = 0;
unprocessed = 1;
// default is pair request
pair = 1;
fix = compute = command = 0;
// default is half neighbor list
half = 1;
full = 0;
full_cluster = 0;
gran = granhistory = 0;
respainner = respamiddle = respaouter = 0;
half_from_full = 0;
// default is every reneighboring
// default is use newton_pair setting in force
// default is encode special bond flags
// default is no auxiliary floating point values
// default is no neighbors of ghosts
// default is no CUDA neighbor list build
// default is no multi-threaded neighbor list build
// default is no Kokkos neighbor list build
+ // default is no Shardlow Splitting Algorithm (SSA) neighbor list build
occasional = 0;
newton = 0;
special = 1;
dnum = 0;
ghost = 0;
cudable = 0;
omp = 0;
intel = 0;
kokkos_host = kokkos_device = 0;
+ ssa = 0;
// default is no copy or skip
copy = 0;
skip = 0;
iskip = NULL;
ijskip = NULL;
otherlist = -1;
}
/* ---------------------------------------------------------------------- */
NeighRequest::~NeighRequest()
{
delete [] iskip;
memory->destroy(ijskip);
}
/* ----------------------------------------------------------------------
archive request params that Neighbor may change after call to identical()
------------------------------------------------------------------------- */
void NeighRequest::archive()
{
half_original = half;
half_from_full_original = half_from_full;
copy_original = copy;
otherlist_original = otherlist;
}
/* ----------------------------------------------------------------------
compare this request to other request
identical means all params set by requestor are the same
compare to original values in other if Neighbor may have changed them
return 1 if identical, 0 if not
------------------------------------------------------------------------- */
int NeighRequest::identical(NeighRequest *other)
{
int same = 1;
// set same = 0 if old list was never processed
// use of requestor_instance and instance counter
// prevents an old fix from being unfix/refix in same memory location
// and appearing to be old, when it is really new
// only needed for classes with persistent neigh lists: Fix, Compute, Pair
if (other->unprocessed) same = 0;
if (requestor != other->requestor) same = 0;
if (requestor_instance != other->requestor_instance) same = 0;
if (id != other->id) same = 0;
if (pair != other->pair) same = 0;
if (fix != other->fix) same = 0;
if (compute != other->compute) same = 0;
if (command != other->command) same = 0;
if (half != other->half_original) same = 0;
if (full != other->full) same = 0;
if (gran != other->gran) same = 0;
if (granhistory != other->granhistory) same = 0;
if (respainner != other->respainner) same = 0;
if (respamiddle != other->respamiddle) same = 0;
if (respaouter != other->respaouter) same = 0;
if (half_from_full != other->half_from_full_original) same = 0;
if (newton != other->newton) same = 0;
if (occasional != other->occasional) same = 0;
if (special != other->special) same = 0;
if (dnum != other->dnum) same = 0;
if (ghost != other->ghost) same = 0;
if (cudable != other->cudable) same = 0;
if (omp != other->omp) same = 0;
if (intel != other->intel) same = 0;
+ if (ssa != other->ssa) same = 0;
if (copy != other->copy_original) same = 0;
if (same_skip(other) == 0) same = 0;
if (otherlist != other->otherlist_original) same = 0;
return same;
}
/* ----------------------------------------------------------------------
compare kind of this request to other request
return 1 if same, 0 if different
------------------------------------------------------------------------- */
int NeighRequest::same_kind(NeighRequest *other)
{
int same = 1;
if (half != other->half) same = 0;
if (full != other->full) same = 0;
if (gran != other->gran) same = 0;
if (granhistory != other->granhistory) same = 0;
if (respainner != other->respainner) same = 0;
if (respamiddle != other->respamiddle) same = 0;
if (respaouter != other->respaouter) same = 0;
if (half_from_full != other->half_from_full) same = 0;
if (newton != other->newton) same = 0;
if (ghost != other->ghost) same = 0;
if (cudable != other->cudable) same = 0;
if (omp != other->omp) same = 0;
if (intel != other->intel) same = 0;
+ if (ssa != other->ssa) same = 0;
return same;
}
/* ----------------------------------------------------------------------
compare skip attributes of this request to other request
return 1 if same, 0 if different
------------------------------------------------------------------------- */
int NeighRequest::same_skip(NeighRequest *other)
{
int i,j;
int same = 1;
if (skip != other->skip) same = 0;
if (skip && other->skip) {
int ntypes = atom->ntypes;
for (i = 1; i <= ntypes; i++)
if (iskip[i] != other->iskip[i]) same = 0;
for (i = 1; i <= ntypes; i++)
for (j = 1; j <= ntypes; j++)
if (ijskip[i][j] != other->ijskip[i][j]) same = 0;
}
return same;
}
/* ----------------------------------------------------------------------
set kind and other values of this request to that of other request
------------------------------------------------------------------------- */
void NeighRequest::copy_request(NeighRequest *other)
{
half = 0;
if (other->half) half = 1;
if (other->full) full = 1;
if (other->gran) gran = 1;
if (other->granhistory) granhistory = 1;
if (other->respainner) respainner = 1;
if (other->respamiddle) respamiddle = 1;
if (other->respaouter) respaouter = 1;
if (other->half_from_full) half_from_full = 1;
newton = other->newton;
dnum = other->dnum;
ghost = other->ghost;
cudable = other->cudable;
omp = other->omp;
intel = other->intel;
+ ssa = other->ssa;
}
diff --git a/src/neigh_request.h b/src/neigh_request.h
index e83523063..047f3d0bd 100644
--- a/src/neigh_request.h
+++ b/src/neigh_request.h
@@ -1,125 +1,129 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifndef LMP_NEIGH_REQUEST_H
#define LMP_NEIGH_REQUEST_H
#include "pointers.h"
namespace LAMMPS_NS {
class NeighRequest : protected Pointers {
public:
void *requestor; // class that made request
int requestor_instance; // instance of that class (only Fix, Compute, Pair)
int id; // ID of request as stored by requestor
// used to track multiple requests from one class
int unprocessed; // 1 when first requested
// 0 after processed by Neighbor class
// which class is requesting the list, one flag is 1, others are 0
int pair; // set by default
int fix;
int compute;
int command;
// kind of list requested, one flag is 1, others are 0
// set by requesting class
int half; // 1 if half neigh list (set by default)
int full; // 1 if full neigh list
int full_cluster; // only used by Kokkos pair styles
int gran; // 1 if granular list
int granhistory; // 1 if granular history list
int respainner; // 1 if a rRESPA inner list
int respamiddle; // 1 if a rRESPA middle list
int respaouter; // 1 if a rRESPA outer list
int half_from_full; // 1 if half list computed from previous full list
// 0 if needed every reneighboring during run
// 1 if occasionally needed by a fix, compute, etc
// set by requesting class
int occasional;
// 0 if use force::newton_pair setting
// 1 if override with pair newton on
// 2 if override with pair newton off
int newton;
// 0 if user of list wants no encoding of special bond flags and all neighs
// 1 if user of list wants special bond flags encoded, set by default
int special;
// number of auxiliary floating point values to store, 0 if none
// set by requesting class
int dnum;
// 1 if also need neighbors of ghosts
int ghost;
// 1 if neighbor list build will be done on GPU
int cudable;
// 1 if using multi-threaded neighbor list build for USER-OMP or USER-INTEL
int omp;
int intel;
// 1 if using Kokkos neighbor build
int kokkos_host;
int kokkos_device;
+ // 1 if using Shardlow Splitting Algorithm (SSA) neighbor list build
+
+ int ssa;
+
// set by neighbor and pair_hybrid after all requests are made
// these settings do not change kind value
int copy; // 1 if this list copied from another list
int skip; // 1 if this list skips atom types from another list
int *iskip; // iskip[i] if atoms of type I are not in list
int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list
int otherlist; // index of other list to copy or skip from
// original params by requestor
// stored to compare against in identical() in case Neighbor changes them
int half_original;
int half_from_full_original;
int copy_original;
int otherlist_original;
// methods
NeighRequest(class LAMMPS *);
~NeighRequest();
void archive();
int identical(NeighRequest *);
int same_kind(NeighRequest *);
int same_skip(NeighRequest *);
void copy_request(NeighRequest *);
};
}
#endif
diff --git a/src/set.cpp b/src/set.cpp
index 182fe528a..a0bd6c2b6 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -1,1054 +1,1066 @@
/* ----------------------------------------------------------------------
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 <math.h>
#include <stdlib.h>
#include <string.h>
#include "set.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "atom_vec_body.h"
#include "domain.h"
#include "region.h"
#include "group.h"
#include "comm.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "input.h"
#include "variable.h"
#include "random_park.h"
#include "math_extra.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,THETA_RANDOM,ANGMOM,OMEGA,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
- MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,INAME,DNAME};
+ MESO_E,MESO_CV,MESO_RHO,SMD_MASS_DENSITY,SMD_CONTACT_RADIUS,DPDTHETA,
+ INAME,DNAME};
#define BIG INT_MAX
/* ---------------------------------------------------------------------- */
void Set::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Set command before simulation box is defined");
if (atom->natoms == 0)
error->all(FLERR,"Set command with no atoms existing");
if (narg < 3) error->all(FLERR,"Illegal set command");
// style and ID info
if (strcmp(arg[0],"atom") == 0) style = ATOM_SELECT;
else if (strcmp(arg[0],"mol") == 0) style = MOL_SELECT;
else if (strcmp(arg[0],"type") == 0) style = TYPE_SELECT;
else if (strcmp(arg[0],"group") == 0) style = GROUP_SELECT;
else if (strcmp(arg[0],"region") == 0) style = REGION_SELECT;
else error->all(FLERR,"Illegal set command");
int n = strlen(arg[1]) + 1;
id = new char[n];
strcpy(id,arg[1]);
select = NULL;
selection(atom->nlocal);
// loop over keyword/value pairs
// call appropriate routine to reset attributes
if (comm->me == 0 && screen) fprintf(screen,"Setting atom values ...\n");
int allcount,origarg;
int iarg = 2;
while (iarg < narg) {
varflag = varflag1 = varflag2 = varflag3 = varflag4 = 0;
count = 0;
origarg = iarg;
if (strcmp(arg[iarg],"type") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else ivalue = force->inumeric(FLERR,arg[iarg+1]);
set(TYPE);
iarg += 2;
} else if (strcmp(arg[iarg],"type/fraction") == 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_FRACTION);
iarg += 4;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (!atom->molecule_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(MOLECULE);
iarg += 2;
} else if (strcmp(arg[iarg],"x") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
set(X);
iarg += 2;
} else if (strcmp(arg[iarg],"y") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
set(Y);
iarg += 2;
} else if (strcmp(arg[iarg],"z") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
set(Z);
iarg += 2;
} else if (strcmp(arg[iarg],"charge") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->q_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(CHARGE);
iarg += 2;
} else if (strcmp(arg[iarg],"mass") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->rmass_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(MASS);
iarg += 2;
} else if (strcmp(arg[iarg],"shape") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else xvalue = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (!atom->ellipsoid_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(SHAPE);
iarg += 4;
} else if (strcmp(arg[iarg],"length") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(LENGTH);
iarg += 2;
} else if (strcmp(arg[iarg],"tri") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->tri_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(TRI);
iarg += 2;
} else if (strcmp(arg[iarg],"dipole") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else xvalue = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (!atom->mu_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(DIPOLE);
iarg += 4;
} else if (strcmp(arg[iarg],"dipole/random") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
dvalue = force->numeric(FLERR,arg[iarg+2]);
if (!atom->mu_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
if (dvalue <= 0.0)
error->all(FLERR,"Invalid dipole length in set command");
setrandom(DIPOLE_RANDOM);
iarg += 3;
} else if (strcmp(arg[iarg],"quat") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else xvalue = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (strstr(arg[iarg+4],"v_") == arg[iarg+4]) varparse(arg[iarg+4],4);
else wvalue = force->numeric(FLERR,arg[iarg+4]);
if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(QUAT);
iarg += 5;
} else if (strcmp(arg[iarg],"quat/random") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
setrandom(QUAT_RANDOM);
iarg += 2;
} else if (strcmp(arg[iarg],"theta") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else {
dvalue = force->numeric(FLERR,arg[iarg+1]);
dvalue *= MY_PI/180.0;
}
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(THETA);
iarg += 2;
} else if (strcmp(arg[iarg],"theta/random") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
set(THETA_RANDOM);
iarg += 2;
} else if (strcmp(arg[iarg],"angmom") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else xvalue = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (!atom->angmom_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(ANGMOM);
iarg += 4;
} else if (strcmp(arg[iarg],"omega") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else xvalue = force->numeric(FLERR,arg[iarg+1]);
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2);
else yvalue = force->numeric(FLERR,arg[iarg+2]);
if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3);
else zvalue = force->numeric(FLERR,arg[iarg+3]);
if (!atom->omega_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(OMEGA);
iarg += 4;
} else if (strcmp(arg[iarg],"diameter") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->radius_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(DIAMETER);
iarg += 2;
} else if (strcmp(arg[iarg],"density") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->rmass_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (dvalue <= 0.0) error->all(FLERR,"Invalid density in set command");
set(DENSITY);
iarg += 2;
} else if (strcmp(arg[iarg],"volume") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->vfrac_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (dvalue <= 0.0) error->all(FLERR,"Invalid volume in set command");
set(VOLUME);
iarg += 2;
} else if (strcmp(arg[iarg],"image") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
ximageflag = yimageflag = zimageflag = 0;
if (strcmp(arg[iarg+1],"NULL") != 0) {
ximageflag = 1;
ximage = force->inumeric(FLERR,arg[iarg+1]);
}
if (strcmp(arg[iarg+2],"NULL") != 0) {
yimageflag = 1;
yimage = force->inumeric(FLERR,arg[iarg+2]);
}
if (strcmp(arg[iarg+3],"NULL") != 0) {
zimageflag = 1;
zimage = force->inumeric(FLERR,arg[iarg+3]);
}
if (ximageflag && ximage && !domain->xperiodic)
error->all(FLERR,
"Cannot set non-zero image flag for non-periodic dimension");
if (yimageflag && yimage && !domain->yperiodic)
error->all(FLERR,
"Cannot set non-zero image flag for non-periodic dimension");
if (zimageflag && zimage && !domain->zperiodic)
error->all(FLERR,
"Cannot set non-zero image flag for non-periodic dimension");
set(IMAGE);
iarg += 4;
} else if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0 || ivalue > atom->nbondtypes)
error->all(FLERR,"Invalid value in set command");
topology(BOND);
iarg += 2;
} else if (strcmp(arg[iarg],"angle") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0 || ivalue > atom->nangletypes)
error->all(FLERR,"Invalid value in set command");
topology(ANGLE);
iarg += 2;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0 || ivalue > atom->ndihedraltypes)
error->all(FLERR,"Invalid value in set command");
topology(DIHEDRAL);
iarg += 2;
} else if (strcmp(arg[iarg],"improper") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
ivalue = force->inumeric(FLERR,arg[iarg+1]);
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0 || ivalue > atom->nimpropertypes)
error->all(FLERR,"Invalid value in set command");
topology(IMPROPER);
iarg += 2;
- } else if (strcmp(arg[iarg],"meso_e") == 0) {
+ } else if (strcmp(arg[iarg],"meso/e") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->e_flag)
- error->all(FLERR,"Cannot set this attribute for this atom style");
+ error->all(FLERR,"Cannot set meso/e for this atom style");
set(MESO_E);
iarg += 2;
- } else if (strcmp(arg[iarg],"meso_cv") == 0) {
+ } else if (strcmp(arg[iarg],"meso/cv") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->cv_flag)
- error->all(FLERR,"Cannot set this attribute for this atom style");
+ error->all(FLERR,"Cannot set meso/cv for this atom style");
set(MESO_CV);
iarg += 2;
- } else if (strcmp(arg[iarg],"meso_rho") == 0) {
+ } else if (strcmp(arg[iarg],"meso/rho") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->rho_flag)
- error->all(FLERR,"Cannot set meso_rho for this atom style");
+ error->all(FLERR,"Cannot set meso/rho for this atom style");
set(MESO_RHO);
iarg += 2;
- } else if (strcmp(arg[iarg],"smd_mass_density") == 0) {
+ } else if (strcmp(arg[iarg],"smd/mass/density") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->smd_flag)
- error->all(FLERR,"Cannot set smd_mass_density for this atom style");
+ error->all(FLERR,"Cannot set smd/mass/density for this atom style");
set(SMD_MASS_DENSITY);
iarg += 2;
- } else if (strcmp(arg[iarg],"smd_contact_radius") == 0) {
+ } else if (strcmp(arg[iarg],"smd/contact/radius") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
if (!atom->smd_flag)
- error->all(FLERR,"Cannot set smd_contact_radius for this atom style");
+ error->all(FLERR,"Cannot set smd/contact/radius "
+ "for this atom style");
set(SMD_CONTACT_RADIUS);
iarg += 2;
+ } else if (strcmp(arg[iarg],"dpd/theta") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+ if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
+ else dvalue = force->numeric(FLERR,arg[iarg+1]);
+ if (!atom->dpd_flag)
+ error->all(FLERR,"Cannot set dpd/theta for this atom style");
+ set(DPDTHETA);
+ iarg += 2;
+
} else if (strstr(arg[iarg],"i_") == arg[iarg]) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else ivalue = force->inumeric(FLERR,arg[iarg+1]);
int flag;
index_custom = atom->find_custom(&arg[iarg][2],flag);
if (index_custom < 0 || flag != 0)
error->all(FLERR,"Set command integer vector does not exist");
set(INAME);
iarg += 2;
} else if (strstr(arg[iarg],"d_") == arg[iarg]) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1);
else dvalue = force->numeric(FLERR,arg[iarg+1]);
int flag;
index_custom = atom->find_custom(&arg[iarg][2],flag);
if (index_custom < 0 || flag != 1)
error->all(FLERR,"Set command floating point vector does not exist");
set(DNAME);
iarg += 2;
} else error->all(FLERR,"Illegal set command");
// statistics
MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world);
if (comm->me == 0) {
if (screen) fprintf(screen," %d settings made for %s\n",
allcount,arg[origarg]);
if (logfile) fprintf(logfile," %d settings made for %s\n",
allcount,arg[origarg]);
}
}
// free local memory
delete [] id;
delete [] select;
}
/* ----------------------------------------------------------------------
select atoms according to ATOM, MOLECULE, TYPE, GROUP, REGION style
n = nlocal or nlocal+nghost depending on keyword
------------------------------------------------------------------------- */
void Set::selection(int n)
{
delete [] select;
select = new int[n];
int nlo,nhi;
if (style == ATOM_SELECT) {
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use set atom with no atom IDs defined");
bigint nlobig,nhibig;
force->boundsbig(id,MAXTAGINT,nlobig,nhibig);
tagint *tag = atom->tag;
for (int i = 0; i < n; i++)
if (tag[i] >= nlobig && tag[i] <= nhibig) select[i] = 1;
else select[i] = 0;
} else if (style == MOL_SELECT) {
if (atom->molecule_flag == 0)
error->all(FLERR,"Cannot use set mol with no molecule IDs defined");
bigint nlobig,nhibig;
force->boundsbig(id,MAXTAGINT,nlobig,nhibig);
tagint *molecule = atom->molecule;
for (int i = 0; i < n; i++)
if (molecule[i] >= nlobig && molecule[i] <= nhibig) select[i] = 1;
else select[i] = 0;
} else if (style == TYPE_SELECT) {
force->bounds(id,atom->ntypes,nlo,nhi);
int *type = atom->type;
for (int i = 0; i < n; i++)
if (type[i] >= nlo && type[i] <= nhi) select[i] = 1;
else select[i] = 0;
} else if (style == GROUP_SELECT) {
int igroup = group->find(id);
if (igroup == -1) error->all(FLERR,"Could not find set group ID");
int groupbit = group->bitmask[igroup];
int *mask = atom->mask;
for (int i = 0; i < n; i++)
if (mask[i] & groupbit) select[i] = 1;
else select[i] = 0;
} else if (style == REGION_SELECT) {
int iregion = domain->find_region(id);
if (iregion == -1) error->all(FLERR,"Set region ID does not exist");
domain->regions[iregion]->prematch();
double **x = atom->x;
for (int i = 0; i < n; i++)
if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
select[i] = 1;
else select[i] = 0;
}
}
/* ----------------------------------------------------------------------
set owned atom properties directly
either scalar or per-atom values from atom-style variable(s)
------------------------------------------------------------------------- */
void Set::set(int keyword)
{
// evaluate atom-style variable(s) if necessary
vec1 = vec2 = vec3 = vec4 = NULL;
if (varflag) {
int nlocal = atom->nlocal;
if (varflag1) {
memory->create(vec1,nlocal,"set:vec1");
input->variable->compute_atom(ivar1,0,vec1,1,0);
}
if (varflag2) {
memory->create(vec2,nlocal,"set:vec2");
input->variable->compute_atom(ivar2,0,vec2,1,0);
}
if (varflag3) {
memory->create(vec3,nlocal,"set:vec3");
input->variable->compute_atom(ivar3,0,vec3,1,0);
}
if (varflag4) {
memory->create(vec4,nlocal,"set:vec4");
input->variable->compute_atom(ivar4,0,vec4,1,0);
}
}
// loop over selected atoms
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body");
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (!select[i]) continue;
// overwrite dvalue, ivalue, xyzw value if variables defined
// else the input script scalar value remains in place
if (varflag) {
if (varflag1) {
dvalue = xvalue = vec1[i];
ivalue = static_cast<int> (dvalue);
}
if (varflag2) yvalue = vec2[i];
if (varflag3) zvalue = vec3[i];
if (varflag4) wvalue = vec4[i];
}
// set values in per-atom arrays
// error check here in case atom-style variables generated bogus value
if (keyword == TYPE) {
if (ivalue <= 0 || ivalue > atom->ntypes)
error->one(FLERR,"Invalid value in set command");
atom->type[i] = ivalue;
}
else if (keyword == MOLECULE) atom->molecule[i] = ivalue;
else if (keyword == X) atom->x[i][0] = dvalue;
else if (keyword == Y) atom->x[i][1] = dvalue;
else if (keyword == Z) atom->x[i][2] = dvalue;
else if (keyword == CHARGE) atom->q[i] = dvalue;
else if (keyword == MASS) {
if (dvalue <= 0.0) error->one(FLERR,"Invalid mass in set command");
atom->rmass[i] = dvalue;
}
else if (keyword == DIAMETER) {
if (dvalue < 0.0) error->one(FLERR,"Invalid diameter in set command");
atom->radius[i] = 0.5 * dvalue;
}
else if (keyword == VOLUME) {
if (dvalue <= 0.0) error->one(FLERR,"Invalid volume in set command");
atom->vfrac[i] = dvalue;
}
else if (keyword == MESO_E) atom->e[i] = dvalue;
else if (keyword == MESO_CV) atom->cv[i] = dvalue;
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
else if (keyword == SMD_MASS_DENSITY) {
// set mass from volume and supplied mass density
atom->rmass[i] = atom->vfrac[i] * dvalue;
}
else if (keyword == SMD_CONTACT_RADIUS) atom->contact_radius[i] = dvalue;
+ else if (keyword == DPDTHETA) atom->dpdTheta[i] = dvalue;
// set shape of ellipsoidal particle
else if (keyword == SHAPE) {
if (xvalue < 0.0 || yvalue < 0.0 || zvalue < 0.0)
error->one(FLERR,"Invalid shape in set command");
if (xvalue > 0.0 || yvalue > 0.0 || zvalue > 0.0) {
if (xvalue == 0.0 || yvalue == 0.0 || zvalue == 0.0)
error->one(FLERR,"Invalid shape in set command");
}
avec_ellipsoid->set_shape(i,0.5*xvalue,0.5*yvalue,0.5*zvalue);
}
// set length of line particle
else if (keyword == LENGTH) {
if (dvalue < 0.0) error->one(FLERR,"Invalid length in set command");
avec_line->set_length(i,dvalue);
}
// set corners of tri particle
else if (keyword == TRI) {
if (dvalue < 0.0) error->one(FLERR,"Invalid length in set command");
avec_tri->set_equilateral(i,dvalue);
}
// set rmass via density
// if radius > 0.0, treat as sphere
// if shape > 0.0, treat as ellipsoid
// if length > 0.0, treat as line
// if area > 0.0, treat as tri
// else set rmass to density directly
else if (keyword == DENSITY) {
if (dvalue <= 0.0) error->one(FLERR,"Invalid density in set command");
if (atom->radius_flag && atom->radius[i] > 0.0)
atom->rmass[i] = 4.0*MY_PI/3.0 *
atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) {
double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape;
atom->rmass[i] = 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue;
} else if (atom->line_flag && atom->line[i] >= 0) {
double length = avec_line->bonus[atom->line[i]].length;
atom->rmass[i] = length * dvalue;
} else if (atom->tri_flag && atom->tri[i] >= 0) {
double *c1 = avec_tri->bonus[atom->tri[i]].c1;
double *c2 = avec_tri->bonus[atom->tri[i]].c2;
double *c3 = avec_tri->bonus[atom->tri[i]].c3;
double c2mc1[3],c3mc1[3];
MathExtra::sub3(c2,c1,c2mc1);
MathExtra::sub3(c3,c1,c3mc1);
double norm[3];
MathExtra::cross3(c2mc1,c3mc1,norm);
double area = 0.5 * MathExtra::len3(norm);
atom->rmass[i] = area * dvalue;
} else atom->rmass[i] = dvalue;
}
// set dipole moment
else if (keyword == DIPOLE) {
double **mu = atom->mu;
mu[i][0] = xvalue;
mu[i][1] = yvalue;
mu[i][2] = zvalue;
mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] +
mu[i][2]*mu[i][2]);
}
// set quaternion orientation of ellipsoid or tri or body particle
// enforce quat rotation vector in z dir for 2d systems
else if (keyword == QUAT) {
double *quat;
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_tri && atom->tri[i] >= 0)
quat = avec_tri->bonus[atom->tri[i]].quat;
else if (avec_body && atom->body[i] >= 0)
quat = avec_body->bonus[atom->body[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
if (domain->dimension == 2 && (xvalue != 0.0 || yvalue != 0.0))
error->one(FLERR,"Cannot set quaternion with xy components "
"for 2d system");
double theta2 = MY_PI2 * wvalue/180.0;
double sintheta2 = sin(theta2);
quat[0] = cos(theta2);
quat[1] = xvalue * sintheta2;
quat[2] = yvalue * sintheta2;
quat[3] = zvalue * sintheta2;
MathExtra::qnormalize(quat);
}
// set theta of line particle
else if (keyword == THETA) {
if (atom->line[i] < 0)
error->one(FLERR,"Cannot set theta for atom that is not a line");
avec_line->bonus[atom->line[i]].theta = dvalue;
}
// set angmom or omega of particle
else if (keyword == ANGMOM) {
atom->angmom[i][0] = xvalue;
atom->angmom[i][1] = yvalue;
atom->angmom[i][2] = zvalue;
}
else if (keyword == OMEGA) {
atom->omega[i][0] = xvalue;
atom->omega[i][1] = yvalue;
atom->omega[i][2] = zvalue;
}
// reset any or all of 3 image flags
else if (keyword == IMAGE) {
int xbox = (atom->image[i] & IMGMASK) - IMGMAX;
int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX;
if (ximageflag) xbox = ximage;
if (yimageflag) ybox = yimage;
if (zimageflag) zbox = zimage;
atom->image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
(((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
}
// set value for custom integer or double vector
else if (keyword == INAME) {
atom->ivector[index_custom][i] = ivalue;
}
else if (keyword == DNAME) {
atom->dvector[index_custom][i] = dvalue;
}
count++;
}
// clear up per-atom memory if allocated
memory->destroy(vec1);
memory->destroy(vec2);
memory->destroy(vec3);
memory->destroy(vec4);
}
/* ----------------------------------------------------------------------
set an owned atom property randomly
set seed based on atom coordinates
make atom result independent of what proc owns it
------------------------------------------------------------------------- */
void Set::setrandom(int keyword)
{
int i;
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body");
RanPark *random = new RanPark(lmp,1);
double **x = atom->x;
int seed = ivalue;
// set fraction of atom types to newtype
if (keyword == TYPE_FRACTION) {
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (select[i]) {
random->reset(seed,x[i]);
if (random->uniform() > fraction) continue;
atom->type[i] = newtype;
count++;
}
// set dipole moments to random orientations in 3d or 2d
// dipole length is determined by dipole type array
} else if (keyword == DIPOLE_RANDOM) {
double **mu = atom->mu;
int nlocal = atom->nlocal;
double msq,scale;
if (domain->dimension == 3) {
for (i = 0; i < nlocal; i++)
if (select[i]) {
random->reset(seed,x[i]);
mu[i][0] = random->uniform() - 0.5;
mu[i][1] = random->uniform() - 0.5;
mu[i][2] = random->uniform() - 0.5;
msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
scale = dvalue/sqrt(msq);
mu[i][0] *= scale;
mu[i][1] *= scale;
mu[i][2] *= scale;
mu[i][3] = dvalue;
count++;
}
} else {
for (i = 0; i < nlocal; i++)
if (select[i]) {
random->reset(seed,x[i]);
mu[i][0] = random->uniform() - 0.5;
mu[i][1] = random->uniform() - 0.5;
mu[i][2] = 0.0;
msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1];
scale = dvalue/sqrt(msq);
mu[i][0] *= scale;
mu[i][1] *= scale;
mu[i][3] = dvalue;
count++;
}
}
// set quaternions to random orientations in 3d and 2d
} else if (keyword == QUAT_RANDOM) {
int nlocal = atom->nlocal;
double *quat;
if (domain->dimension == 3) {
double s,t1,t2,theta1,theta2;
for (i = 0; i < nlocal; i++)
if (select[i]) {
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_tri && atom->tri[i] >= 0)
quat = avec_tri->bonus[atom->tri[i]].quat;
else if (avec_body && atom->body[i] >= 0)
quat = avec_body->bonus[atom->body[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
random->reset(seed,x[i]);
s = random->uniform();
t1 = sqrt(1.0-s);
t2 = sqrt(s);
theta1 = 2.0*MY_PI*random->uniform();
theta2 = 2.0*MY_PI*random->uniform();
quat[0] = cos(theta2)*t2;
quat[1] = sin(theta1)*t1;
quat[2] = cos(theta1)*t1;
quat[3] = sin(theta2)*t2;
count++;
}
} else {
double theta2;
for (i = 0; i < nlocal; i++)
if (select[i]) {
if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
else if (avec_body && atom->body[i] >= 0)
quat = avec_body->bonus[atom->body[i]].quat;
else
error->one(FLERR,"Cannot set quaternion for atom that has none");
random->reset(seed,x[i]);
theta2 = MY_PI*random->uniform();
quat[0] = cos(theta2);
quat[1] = 0.0;
quat[2] = 0.0;
quat[3] = sin(theta2);
count++;
}
}
// set theta to random orientation in 2d
} else if (keyword == THETA_RANDOM) {
int nlocal = atom->nlocal;
double theta;
for (i = 0; i < nlocal; i++) {
if (select[i]) {
if (atom->line[i] < 0)
error->one(FLERR,"Cannot set theta for atom that is not a line");
random->reset(seed,x[i]);
avec_line->bonus[atom->line[i]].theta = MY_2PI*random->uniform();
count++;
}
}
}
delete random;
}
/* ---------------------------------------------------------------------- */
void Set::topology(int keyword)
{
int m,atom1,atom2,atom3,atom4;
// error check
if (atom->molecular == 2)
error->all(FLERR,"Cannot set bond topology types for atom style template");
// border swap to acquire ghost atom info
// enforce PBC before in case atoms are outside box
// init entire system since comm->exchange is done
// comm::init needs neighbor::init needs pair::init needs kspace::init, etc
if (comm->me == 0 && screen) fprintf(screen," system init for set ...\n");
lmp->init();
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);
// select both owned and ghost atoms
selection(atom->nlocal + atom->nghost);
// for BOND, each of 2 atoms must be in group
if (keyword == BOND) {
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_bond[i]; m++) {
atom1 = atom->map(atom->bond_atom[i][m]);
if (atom1 == -1) error->one(FLERR,"Bond atom missing in set command");
if (select[i] && select[atom1]) {
atom->bond_type[i][m] = ivalue;
count++;
}
}
}
// for ANGLE, each of 3 atoms must be in group
if (keyword == ANGLE) {
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_angle[i]; m++) {
atom1 = atom->map(atom->angle_atom1[i][m]);
atom2 = atom->map(atom->angle_atom2[i][m]);
atom3 = atom->map(atom->angle_atom3[i][m]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1)
error->one(FLERR,"Angle atom missing in set command");
if (select[atom1] && select[atom2] && select[atom3]) {
atom->angle_type[i][m] = ivalue;
count++;
}
}
}
// for DIHEDRAL, each of 4 atoms must be in group
if (keyword == DIHEDRAL) {
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_dihedral[i]; m++) {
atom1 = atom->map(atom->dihedral_atom1[i][m]);
atom2 = atom->map(atom->dihedral_atom2[i][m]);
atom3 = atom->map(atom->dihedral_atom3[i][m]);
atom4 = atom->map(atom->dihedral_atom4[i][m]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
error->one(FLERR,"Dihedral atom missing in set command");
if (select[atom1] && select[atom2] && select[atom3] && select[atom4]) {
atom->dihedral_type[i][m] = ivalue;
count++;
}
}
}
// for IMPROPER, each of 4 atoms must be in group
if (keyword == IMPROPER) {
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
for (m = 0; m < atom->num_improper[i]; m++) {
atom1 = atom->map(atom->improper_atom1[i][m]);
atom2 = atom->map(atom->improper_atom2[i][m]);
atom3 = atom->map(atom->improper_atom3[i][m]);
atom4 = atom->map(atom->improper_atom4[i][m]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
error->one(FLERR,"Improper atom missing in set command");
if (select[atom1] && select[atom2] && select[atom3] && select[atom4]) {
atom->improper_type[i][m] = ivalue;
count++;
}
}
}
}
/* ---------------------------------------------------------------------- */
void Set::varparse(char *name, int m)
{
varflag = 1;
name = &name[2];
int n = strlen(name) + 1;
char *str = new char[n];
strcpy(str,name);
int ivar = input->variable->find(str);
delete [] str;
if (ivar < 0)
error->all(FLERR,"Variable name for set command does not exist");
if (!input->variable->atomstyle(ivar))
error->all(FLERR,"Variable for set command is invalid style");
if (m == 1) {
varflag1 = 1; ivar1 = ivar;
} else if (m == 2) {
varflag2 = 1; ivar2 = ivar;
} else if (m == 3) {
varflag3 = 1; ivar3 = ivar;
} else if (m == 4) {
varflag4 = 1; ivar4 = ivar;
}
}