ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
4c0b4154879e49f43fc2926d3d62369c732d8653
lammps/src/atom.cpp
Axel Kohlmeyer 248b7383c7 Merge remote-tracking branch 'github/develop' into collected-small-changes 
# Conflicts:
#	src/DIPOLE/pair_lj_cut_dipole_long.h
#	src/OPENMP/fix_nvt_sllod_omp.h
#	src/PHONON/fix_phonon.h
#	src/SPIN/pair_spin_dipole_cut.h
#	src/SPIN/pair_spin_dipole_long.h
#	src/STUBS/mpi.h
#	src/UEF/fix_nh_uef.h
#	src/fix_nvt_sllod.h
#	tools/coding_standard/homepage.py
2022-10-24 15:50:50 -04:00

3126 lines
104 KiB
C++
Raw Blame History

// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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 "atom.h"
#include "atom_vec.h"
#include "style_atom.h" // IWYU pragma: keep
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "input.h"
#include "label_map.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "molecule.h"
#include "neighbor.h"
#include "tokenizer.h"
#include "update.h"
#include "variable.h"
#include "library.h"
#include <algorithm>
#include <cstring>
#ifdef LMP_GPU
#include "fix_gpu.h"
#include <cmath>
#endif
using namespace LAMMPS_NS;
using namespace MathConst;
#define DELTA 1
#define EPSILON 1.0e-6
#define MAXLINE 256
/* ----------------------------------------------------------------------
one instance per AtomVec style in style_atom.h
------------------------------------------------------------------------- */
template <typename T> static AtomVec *avec_creator(LAMMPS *_lmp)
{
return new T(_lmp);
}
/* ---------------------------------------------------------------------- */
/** \class LAMMPS_NS::Atom
* \brief Class to provide access to atom data
\verbatim embed:rst
The Atom class provides access to atom style related global settings and
per-atom data that is stored with atoms and migrates with them from
sub-domain to sub-domain as atoms move around. This includes topology
data, which is stored with either one specific atom or all atoms involved
depending on the settings of the :doc:`newton command <newton>`.
The actual per-atom data is allocated and managed by one of the various
classes derived from the AtomVec class as determined by
the :doc:`atom_style command <atom_style>`. The pointers in the Atom class
are updated by the AtomVec class as needed.
\endverbatim
*/
/** Atom class constructor
*
* This resets and initialized all kinds of settings,
* parameters, and pointer variables for per-atom arrays.
* This also initializes the factory for creating
* instances of classes derived from the AtomVec base
* class, which correspond to the selected atom style.
*
* \param _lmp pointer to the base LAMMPS class */
Atom::Atom(LAMMPS *_lmp) : Pointers(_lmp)
{
natoms = 0;
nlocal = nghost = nmax = 0;
ntypes = 0;
nellipsoids = nlines = ntris = nbodies = 0;
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
nbonds = nangles = ndihedrals = nimpropers = 0;
firstgroupname = nullptr;
sortfreq = 1000;
nextsort = 0;
userbinsize = 0.0;
maxbin = maxnext = 0;
binhead = nullptr;
next = permute = nullptr;
// --------------------------------------------------------------------
// 1st customization section: customize by adding new per-atom variables
tag = nullptr;
type = mask = nullptr;
image = nullptr;
x = v = f = nullptr;
// charged and dipolar particles
q = nullptr;
mu = nullptr;
// finite-size particles
omega = angmom = torque = nullptr;
radius = rmass = nullptr;
ellipsoid = line = tri = body = nullptr;
quat = nullptr;
// molecular systems
molecule = nullptr;
molindex = molatom = nullptr;
bond_per_atom = extra_bond_per_atom = 0;
num_bond = nullptr;
bond_type = nullptr;
bond_atom = nullptr;
angle_per_atom = extra_angle_per_atom = 0;
num_angle = nullptr;
angle_type = nullptr;
angle_atom1 = angle_atom2 = angle_atom3 = nullptr;
dihedral_per_atom = extra_dihedral_per_atom = 0;
num_dihedral = nullptr;
dihedral_type = nullptr;
dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = nullptr;
improper_per_atom = extra_improper_per_atom = 0;
num_improper = nullptr;
improper_type = nullptr;
improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = nullptr;
maxspecial = 1;
nspecial = nullptr;
special = nullptr;
// PERI package
vfrac = s0 = nullptr;
x0 = nullptr;
// SPIN package
sp = fm = fm_long = nullptr;
// EFF and AWPMD packages
spin = nullptr;
eradius = ervel = erforce = nullptr;
ervelforce = nullptr;
cs = csforce = vforce = nullptr;
etag = nullptr;
// CG-DNA package
id5p = nullptr;
// DPD-REACT package
uCond = uMech = uChem = uCG = uCGnew = nullptr;
duChem = dpdTheta = nullptr;
// MESO package
cc = cc_flux = nullptr;
edpd_temp = edpd_flux = edpd_cv = nullptr;
// MESONT package
length = nullptr;
buckling = nullptr;
bond_nt = nullptr;
// MACHDYN package
contact_radius = nullptr;
smd_data_9 = nullptr;
smd_stress = nullptr;
eff_plastic_strain = nullptr;
eff_plastic_strain_rate = nullptr;
damage = nullptr;
// SPH package
rho = drho = esph = desph = cv = nullptr;
vest = nullptr;
// AMOEBA package
maxspecial15 = 1;
nspecial15 = nullptr;
special15 = nullptr;
// DIELECTRIC package
area = ed = em = epsilon = curvature = q_unscaled = nullptr;
// end of customization section
// --------------------------------------------------------------------
// user-defined molecules
nmolecule = 0;
molecules = nullptr;
// type labels
lmap = nullptr;
types_style = NUMERIC;
// custom atom arrays
nivector = ndvector = niarray = ndarray = 0;
ivector = nullptr;
dvector = nullptr;
iarray = nullptr;
darray = nullptr;
icols = dcols = nullptr;
ivname = dvname = ianame = daname = nullptr;
// initialize atom style and array existence flags
set_atomflag_defaults();
// initialize peratom data structure
peratom_create();
// ntype-length arrays
mass = nullptr;
mass_setflag = nullptr;
// callback lists & extra restart info
nextra_grow = nextra_restart = nextra_border = 0;
extra_grow = extra_restart = extra_border = nullptr;
nextra_grow_max = nextra_restart_max = nextra_border_max = 0;
nextra_store = 0;
extra = nullptr;
// default atom ID and mapping values
tag_enable = 1;
map_style = map_user = MAP_NONE;
map_tag_max = -1;
map_maxarray = map_nhash = map_nbucket = -1;
max_same = 0;
sametag = nullptr;
map_array = nullptr;
map_bucket = nullptr;
map_hash = nullptr;
unique_tags = nullptr;
reset_image_flag[0] = reset_image_flag[1] = reset_image_flag[2] = false;
atom_style = nullptr;
avec = nullptr;
avec_map = new AtomVecCreatorMap();
#define ATOM_CLASS
#define AtomStyle(key,Class) \
(*avec_map)[#key] = &avec_creator<Class>;
#include "style_atom.h" // IWYU pragma: keep
#undef AtomStyle
#undef ATOM_CLASS
}
/* ---------------------------------------------------------------------- */
Atom::~Atom()
{
delete[] atom_style;
delete avec;
delete avec_map;
delete[] firstgroupname;
memory->destroy(binhead);
memory->destroy(next);
memory->destroy(permute);
memory->destroy(tag);
memory->destroy(type);
memory->destroy(mask);
memory->destroy(image);
memory->destroy(x);
memory->destroy(v);
memory->destroy(f);
// delete custom atom arrays
for (int i = 0; i < nivector; i++) {
delete[] ivname[i];
memory->destroy(ivector[i]);
}
for (int i = 0; i < ndvector; i++) {
delete[] dvname[i];
if (dvector) // (needed for Kokkos)
memory->destroy(dvector[i]);
}
for (int i = 0; i < niarray; i++) {
delete[] ianame[i];
memory->destroy(iarray[i]);
}
for (int i = 0; i < ndarray; i++) {
delete[] daname[i];
memory->destroy(darray[i]);
}
memory->sfree(ivname);
memory->sfree(dvname);
memory->sfree(ianame);
memory->sfree(daname);
memory->sfree(ivector);
memory->sfree(dvector);
memory->sfree(iarray);
memory->sfree(darray);
memory->sfree(icols);
memory->sfree(dcols);
// delete user-defined molecules
for (int i = 0; i < nmolecule; i++) delete molecules[i];
memory->sfree(molecules);
// delete label map
delete lmap;
// 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
Atom::map_delete();
delete unique_tags;
}
/* ----------------------------------------------------------------------
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)
firstgroupname = utils::strdup(old->firstgroupname);
}
/* ----------------------------------------------------------------------
one-time creation of peratom data structure
------------------------------------------------------------------------- */
void Atom::peratom_create()
{
peratom.clear();
// --------------------------------------------------------------------
// 2nd customization section: add peratom variables here, order does not matter
// register tagint & imageint variables as INT or BIGINT
int tagintsize = INT;
if (sizeof(tagint) == 8) tagintsize = BIGINT;
int imageintsize = INT;
if (sizeof(imageint) == 8) imageintsize = BIGINT;
add_peratom("id",&tag,tagintsize,0);
add_peratom("type",&type,INT,0);
add_peratom("mask",&mask,INT,0);
add_peratom("image",&image,imageintsize,0);
add_peratom("x",&x,DOUBLE,3);
add_peratom("v",&v,DOUBLE,3);
add_peratom("f",&f,DOUBLE,3,1); // set per-thread flag
add_peratom("rmass",&rmass,DOUBLE,0);
add_peratom("q",&q,DOUBLE,0);
add_peratom("mu",&mu,DOUBLE,4);
add_peratom("mu3",&mu,DOUBLE,3); // just first 3 values of mu[4]
// finite size particles
add_peratom("radius",&radius,DOUBLE,0);
add_peratom("omega",&omega,DOUBLE,3);
add_peratom("torque",&torque,DOUBLE,3,1); // set per-thread flag
add_peratom("angmom",&angmom,DOUBLE,3);
add_peratom("ellipsoid",&ellipsoid,INT,0);
add_peratom("line",&line,INT,0);
add_peratom("tri",&tri,INT,0);
add_peratom("body",&body,INT,0);
// BPM package
add_peratom("quat",&quat,DOUBLE,4);
// MOLECULE package
add_peratom("molecule",&molecule,tagintsize,0);
add_peratom("molindex",&molindex,INT,0);
add_peratom("molatom",&molatom,INT,0);
add_peratom("nspecial",&nspecial,INT,3);
add_peratom_vary("special",&special,tagintsize,&maxspecial,&nspecial,3);
add_peratom("num_bond",&num_bond,INT,0);
add_peratom_vary("bond_type",&bond_type,INT,&bond_per_atom,&num_bond);
add_peratom_vary("bond_atom",&bond_atom,tagintsize,&bond_per_atom,&num_bond);
add_peratom("num_angle",&num_angle,INT,0);
add_peratom_vary("angle_type",&angle_type,INT,&angle_per_atom,&num_angle);
add_peratom_vary("angle_atom1",&angle_atom1,tagintsize,
&angle_per_atom,&num_angle);
add_peratom_vary("angle_atom2",&angle_atom2,tagintsize,
&angle_per_atom,&num_angle);
add_peratom_vary("angle_atom3",&angle_atom3,tagintsize,
&angle_per_atom,&num_angle);
add_peratom("num_dihedral",&num_dihedral,INT,0);
add_peratom_vary("dihedral_type",&dihedral_type,INT,
&dihedral_per_atom,&num_dihedral);
add_peratom_vary("dihedral_atom1",&dihedral_atom1,tagintsize,
&dihedral_per_atom,&num_dihedral);
add_peratom_vary("dihedral_atom2",&dihedral_atom2,tagintsize,
&dihedral_per_atom,&num_dihedral);
add_peratom_vary("dihedral_atom3",&dihedral_atom3,tagintsize,
&dihedral_per_atom,&num_dihedral);
add_peratom_vary("dihedral_atom4",&dihedral_atom4,tagintsize,
&dihedral_per_atom,&num_dihedral);
add_peratom("num_improper",&num_improper,INT,0);
add_peratom_vary("improper_type",&improper_type,INT,
&improper_per_atom,&num_improper);
add_peratom_vary("improper_atom1",&improper_atom1,tagintsize,
&improper_per_atom,&num_improper);
add_peratom_vary("improper_atom2",&improper_atom2,tagintsize,
&improper_per_atom,&num_improper);
add_peratom_vary("improper_atom3",&improper_atom3,tagintsize,
&improper_per_atom,&num_improper);
add_peratom_vary("improper_atom4",&improper_atom4,tagintsize,
&improper_per_atom,&num_improper);
// PERI package
add_peratom("vfrac",&vfrac,DOUBLE,0);
add_peratom("s0",&s0,DOUBLE,0);
add_peratom("x0",&x0,DOUBLE,3);
// SPIN package
add_peratom("sp",&sp,DOUBLE,4);
add_peratom("fm",&fm,DOUBLE,3,1);
add_peratom("fm_long",&fm_long,DOUBLE,3,1);
// EFF package
add_peratom("espin",&spin,INT,0);
add_peratom("eradius",&eradius,DOUBLE,0);
add_peratom("ervel",&ervel,DOUBLE,0);
add_peratom("erforce",&erforce,DOUBLE,0,1); // set per-thread flag
// AWPMD package
add_peratom("cs",&cs,DOUBLE,2);
add_peratom("csforce",&csforce,DOUBLE,2);
add_peratom("vforce",&vforce,DOUBLE,3);
add_peratom("ervelforce",&ervelforce,DOUBLE,0);
add_peratom("etag",&etag,INT,0);
// CG-DNA package
add_peratom("id5p",&id5p,tagintsize,0);
// DPD-REACT package
add_peratom("dpdTheta",&dpdTheta,DOUBLE,0);
add_peratom("uCond",&uCond,DOUBLE,0);
add_peratom("uMech",&uMech,DOUBLE,0);
add_peratom("uChem",&uChem,DOUBLE,0);
add_peratom("uCG",&uCG,DOUBLE,0);
add_peratom("uCGnew",&uCGnew,DOUBLE,0);
add_peratom("duChem",&duChem,DOUBLE,0);
// MESO package
add_peratom("edpd_cv",&edpd_cv,DOUBLE,0);
add_peratom("edpd_temp",&edpd_temp,DOUBLE,0);
add_peratom("vest_temp",&vest_temp,DOUBLE,0);
add_peratom("edpd_flux",&edpd_flux,DOUBLE,0,1); // set per-thread flag
add_peratom("cc",&cc,DOUBLE,1);
add_peratom("cc_flux",&cc_flux,DOUBLE,1,1); // set per-thread flag
// MESONT package
add_peratom("length",&length,DOUBLE,0);
add_peratom("buckling",&buckling,INT,0);
add_peratom("bond_nt",&bond_nt,tagintsize,2);
// SPH package
add_peratom("rho",&rho,DOUBLE,0);
add_peratom("drho",&drho,DOUBLE,0,1); // set per-thread flag
add_peratom("esph",&esph,DOUBLE,0);
add_peratom("desph",&desph,DOUBLE,0,1); // set per-thread flag
add_peratom("vest",&vest,DOUBLE,3);
add_peratom("cv",&cv,DOUBLE,0);
// MACHDYN package
add_peratom("contact_radius",&contact_radius,DOUBLE,0);
add_peratom("smd_data_9",&smd_data_9,DOUBLE,1);
add_peratom("smd_stress",&smd_stress,DOUBLE,1);
add_peratom("eff_plastic_strain",&eff_plastic_strain,DOUBLE,0);
add_peratom("eff_plastic_strain_rate",&eff_plastic_strain_rate,DOUBLE,0);
add_peratom("damage",&damage,DOUBLE,0);
// AMOEBA package
add_peratom("nspecial15",&nspecial15,INT,0);
add_peratom_vary("special15",&special15,tagintsize,&maxspecial15,&nspecial15,0);
// DIELECTRIC package
add_peratom("area",&area,DOUBLE,0);
add_peratom("ed",&ed,DOUBLE,0);
add_peratom("em",&em,DOUBLE,0);
add_peratom("epsilon",&epsilon,DOUBLE,0);
add_peratom("curvature",&curvature,DOUBLE,0);
add_peratom("q_unscaled",&q_unscaled,DOUBLE,0);
// end of customization section
// --------------------------------------------------------------------
}
/* ----------------------------------------------------------------------
add info for a single per-atom vector/array to PerAtom data struct
cols = 0: per-atom vector
cols = N: static per-atom array with N columns
use add_peratom_vary() when column count varies per atom
------------------------------------------------------------------------- */
void Atom::add_peratom(const std::string &name, void *address,
int datatype, int cols, int threadflag)
{
PerAtom item = {name, address, nullptr, nullptr, datatype, cols, 0, threadflag};
peratom.push_back(item);
}
/* ----------------------------------------------------------------------
change the column count of an existing peratom array entry
allows atom_style to specify column count as an argument
see atom_style tdpd as an example
------------------------------------------------------------------------- */
void Atom::add_peratom_change_columns(const std::string &name, int cols)
{
auto match = std::find_if(peratom.begin(), peratom.end(),
[&name] (const PerAtom &p) { return p.name == name; });
if (match != peratom.end()) (*match).cols = cols;
else error->all(FLERR,"Could not find per-atom array name {} for column change", name);
}
/* ----------------------------------------------------------------------
add info for a single per-atom array to PerAtom data struct
cols = address of int variable with max columns per atom
for collength = 0:
length = address of peratom vector with column count per atom
e.g. num_bond
for collength = N:
length = address of peratom array with column count per atom
collength = index of column (1 to N) in peratom array with count
e.g. nspecial
------------------------------------------------------------------------- */
void Atom::add_peratom_vary(const std::string &name, void *address,
int datatype, int *cols, void *length, int collength)
{
PerAtom item = {name, address, length, cols, datatype, -1, collength, 0};
peratom.push_back(item);
}
/* ----------------------------------------------------------------------
add info for a single per-atom array to PerAtom data struct
------------------------------------------------------------------------- */
void Atom::set_atomflag_defaults()
{
// --------------------------------------------------------------------
// 3rd customization section: customize by adding new flag
// identical list as 2nd customization in atom.h
labelmapflag = 0;
sphere_flag = ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
quat_flag = 0;
peri_flag = electron_flag = 0;
wavepacket_flag = sph_flag = 0;
molecule_flag = molindex_flag = molatom_flag = 0;
q_flag = mu_flag = 0;
rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0;
rho_flag = esph_flag = cv_flag = vest_flag = 0;
dpd_flag = edpd_flag = tdpd_flag = 0;
sp_flag = 0;
x0_flag = 0;
smd_flag = damage_flag = 0;
mesont_flag = 0;
contact_radius_flag = smd_data_9_flag = smd_stress_flag = 0;
eff_plastic_strain_flag = eff_plastic_strain_rate_flag = 0;
nspecial15_flag = 0;
pdscale = 1.0;
}
/* ----------------------------------------------------------------------
create an AtomVec style
called from lammps.cpp, input script, restart file, replicate
------------------------------------------------------------------------- */
void Atom::create_avec(const std::string &style, int narg, char **arg, int trysuffix)
{
delete[] atom_style;
if (avec) delete avec;
atom_style = nullptr;
avec = nullptr;
// unset atom style and array existence flags
// may have been set by old avec
set_atomflag_defaults();
// 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) {
std::string estyle = style + "/";
if (sflag == 1) estyle += lmp->suffix;
else estyle += lmp->suffix2;
atom_style = utils::strdup(estyle);
} else {
atom_style = utils::strdup(style);
}
// if molecular system:
// atom IDs must be defined
// force atom map to be created
// map style will be reset to array vs hash to by map_init()
molecular = avec->molecular;
if ((molecular != Atom::ATOMIC) && (tag_enable == 0))
error->all(FLERR,"Atom IDs must be used for molecular systems");
if (molecular != Atom::ATOMIC) map_style = MAP_YES;
}
/* ----------------------------------------------------------------------
generate an AtomVec class, first with suffix appended
------------------------------------------------------------------------- */
AtomVec *Atom::new_avec(const std::string &style, int trysuffix, int &sflag)
{
if (trysuffix && lmp->suffix_enable) {
if (lmp->suffix) {
sflag = 1;
std::string estyle = style + "/" + lmp->suffix;
if (avec_map->find(estyle) != avec_map->end()) {
AtomVecCreator &avec_creator = (*avec_map)[estyle];
return avec_creator(lmp);
}
}
if (lmp->suffix2) {
sflag = 2;
std::string estyle = style + "/" + lmp->suffix2;
if (avec_map->find(estyle) != avec_map->end()) {
AtomVecCreator &avec_creator = (*avec_map)[estyle];
return avec_creator(lmp);
}
}
}
sflag = 0;
if (avec_map->find(style) != avec_map->end()) {
AtomVecCreator &avec_creator = (*avec_map)[style];
return avec_creator(lmp);
}
error->all(FLERR,utils::check_packages_for_style("atom",style,lmp));
return nullptr;
}
/* ---------------------------------------------------------------------- */
void Atom::init()
{
// delete extra array since it doesn't persist past first run
if (nextra_store) {
memory->destroy(extra);
extra = nullptr;
nextra_store = 0;
}
// check arrays that are atom type in length
check_mass(FLERR);
// setup of firstgroup
if (firstgroupname) {
firstgroup = group->find(firstgroupname);
if (firstgroup < 0)
error->all(FLERR,"Could not find atom_modify first group ID {}", firstgroupname);
} 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();
}
/* ---------------------------------------------------------------------- */
std::string Atom::get_style()
{
std::string retval = atom_style;
if (retval == "hybrid") {
auto avec_hybrid = dynamic_cast<AtomVecHybrid *>(avec);
if (avec_hybrid) {
for (int i = 0; i < avec_hybrid->nstyles; i++) {
retval += ' ';
retval += avec_hybrid->keywords[i];
}
}
}
return retval;
}
/* ----------------------------------------------------------------------
return ptr to AtomVec class if matches style or to matching hybrid sub-class
return nullptr 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) {
auto avec_hybrid = dynamic_cast<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 nullptr;
}
/* ----------------------------------------------------------------------
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) utils::missing_cmd_args(FLERR, "atom_modify", error);
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"id") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "atom_modify id", error);
if (domain->box_exist)
error->all(FLERR,"Atom_modify id command after simulation box is defined");
tag_enable = utils::logical(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"map") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "atom_modify map", error);
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 if (strcmp(arg[iarg+1],"yes") == 0) map_user = 3;
else error->all(FLERR,"Illegal atom_modify map command argument {}", arg[iarg+1]);
map_style = map_user;
iarg += 2;
} else if (strcmp(arg[iarg],"first") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "atom_modify first", error);
if (strcmp(arg[iarg+1],"all") == 0) {
delete[] firstgroupname;
firstgroupname = nullptr;
} else {
firstgroupname = utils::strdup(arg[iarg+1]);
sortfreq = 0;
}
iarg += 2;
} else if (strcmp(arg[iarg],"sort") == 0) {
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "atom_modify sort", error);
sortfreq = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
userbinsize = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (sortfreq < 0) error->all(FLERR,"Illegal atom_modify sort frequency {}", sortfreq);
if (userbinsize < 0.0) error->all(FLERR,"Illegal atom_modify sort bin size {}", userbinsize);
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 argument: {}", arg[iarg]);
}
}
/* ----------------------------------------------------------------------
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 are negative");
if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs are too big");
if (maxall > 0 && minall == 0)
error->all(FLERR,"One or more atom IDs are 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(&notag,&notag_total,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (notag_total >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID {}", MAXTAGINT);
bigint notag_sum;
MPI_Scan(&notag,&notag_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;
}
/* ----------------------------------------------------------------------
check that bonus data settings are valid
error if number of atoms with ellipsoid/line/tri/body flags
are consistent with global setting.
------------------------------------------------------------------------- */
void Atom::bonus_check()
{
bigint local_ellipsoids = 0, local_lines = 0, local_tris = 0;
bigint local_bodies = 0, num_global;
for (int i = 0; i < nlocal; ++i) {
if (ellipsoid && (ellipsoid[i] >=0)) ++local_ellipsoids;
if (line && (line[i] >=0)) ++local_lines;
if (tri && (tri[i] >=0)) ++local_tris;
if (body && (body[i] >=0)) ++local_bodies;
}
MPI_Allreduce(&local_ellipsoids,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (nellipsoids != num_global)
error->all(FLERR,"Inconsistent 'ellipsoids' header value and number of "
"atoms with enabled ellipsoid flags");
MPI_Allreduce(&local_lines,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (nlines != num_global)
error->all(FLERR,"Inconsistent 'lines' header value and number of "
"atoms with enabled line flags");
MPI_Allreduce(&local_tris,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (ntris != num_global)
error->all(FLERR,"Inconsistent 'tris' header value and number of "
"atoms with enabled tri flags");
MPI_Allreduce(&local_bodies,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (nbodies != num_global)
error->all(FLERR,"Inconsistent 'bodies' header value and number of "
"atoms with enabled body flags");
}
/* ----------------------------------------------------------------------
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 = nullptr;
atom->bond_atom = nullptr;
memory->destroy(atom->angle_type);
memory->destroy(atom->angle_atom1);
memory->destroy(atom->angle_atom2);
memory->destroy(atom->angle_atom3);
atom->angle_type = nullptr;
atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = nullptr;
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 = nullptr;
atom->dihedral_atom1 = atom->dihedral_atom2 =
atom->dihedral_atom3 = atom->dihedral_atom4 = nullptr;
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 = nullptr;
atom->improper_atom1 = atom->improper_atom2 =
atom->improper_atom3 = atom->improper_atom4 = nullptr;
}
/* ----------------------------------------------------------------------
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, tagint mol_offset,
int type_offset, int shiftflag, double *shift,
int labelflag, int *ilabel)
{
int xptr,iptr;
imageint imagedata;
double xdata[3],lamda[3];
double *coord;
char *next;
std::string typestr;
auto location = "Atoms section of data file";
// use the first line to detect and validate the number of words/tokens per line
next = strchr(buf,'\n');
if (!next) error->all(FLERR, "Missing data in {}", location);
*next = '\0';
auto values = Tokenizer(buf).as_vector();
int nwords = values.size();
for (std::size_t i = 0; i < values.size(); ++i) {
if (utils::strmatch(values[i], "^#")) {
nwords = i;
break;
}
}
if ((nwords != avec->size_data_atom) && (nwords != avec->size_data_atom + 3))
error->all(FLERR,"Incorrect format in {}: {}", location, utils::trim(buf));
*next = '\n';
// 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 != Comm::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');
if (!next) error->all(FLERR, "Missing data in {}", location);
*next = '\0';
auto values = Tokenizer(buf).as_vector();
int nvalues = values.size();
if ((nvalues == 0) || (utils::strmatch(values[0],"^#.*"))) {
// skip over empty or comment lines
} else if ((nvalues < nwords) ||
((nvalues > nwords) && (!utils::strmatch(values[nwords],"^#")))) {
error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
} else {
int imx = 0, imy = 0, imz = 0;
if (imageflag) {
imx = utils::inumeric(FLERR,values[iptr],false,lmp);
imy = utils::inumeric(FLERR,values[iptr+1],false,lmp);
imz = utils::inumeric(FLERR,values[iptr+2],false,lmp);
if ((domain->dimension == 2) && (imz != 0))
error->all(FLERR,"Z-direction image flag must be 0 for 2d-systems");
if ((!domain->xperiodic) && (imx != 0)) { reset_image_flag[0] = true; imx = 0; }
if ((!domain->yperiodic) && (imy != 0)) { reset_image_flag[1] = true; imy = 0; }
if ((!domain->zperiodic) && (imz != 0)) { reset_image_flag[2] = true; imz = 0; }
}
imagedata = ((imageint) (imx + IMGMAX) & IMGMASK) |
(((imageint) (imy + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (imz + IMGMAX) & IMGMASK) << IMG2BITS);
xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp);
xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp);
xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp);
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,typestr);
typestr = utils::utf8_subst(typestr);
if (id_offset) tag[nlocal-1] += id_offset;
if (mol_offset) molecule[nlocal-1] += mol_offset;
switch (utils::is_type(typestr)) {
case 0: { // numeric
int itype = utils::inumeric(FLERR, typestr, true, lmp) + type_offset;
if ((itype < 1) || (itype > ntypes))
error->one(FLERR, "Invalid atom type {} in {}: {}", itype, location,
utils::trim(buf));
type[nlocal - 1] = itype;
if (labelflag) type[nlocal - 1] = ilabel[itype - 1];
break;
}
case 1: { // type label
if (!labelmapflag)
error->one(FLERR, "Invalid line in {}: {}", location, utils::trim(buf));
type[nlocal - 1] = lmap->find(typestr, Atom::ATOM);
if (type[nlocal - 1] == -1)
error->one(FLERR, "Invalid line in {}: {}", location, utils::trim(buf));
break;
}
default: // invalid
error->one(FLERR, "Invalid line in {}: {}", location, utils::trim(buf));
break;
}
if (type[nlocal-1] <= 0 || type[nlocal-1] > ntypes)
error->one(FLERR,"Invalid atom type {} in {}", location, typestr);
}
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
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 m;
char *next;
// 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');
if (!next) error->all(FLERR, "Missing data in Velocities section of data file");
*next = '\0';
auto values = Tokenizer(utils::trim_comment(buf)).as_vector();
if (values.size() == 0) {
// skip over empty or comment lines
} else if ((int)values.size() != avec->size_data_vel) {
error->all(FLERR, "Incorrect velocity format in data file: {}", utils::trim(buf));
} else {
tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset;
if (tagdata <= 0 || tagdata > map_tag_max)
error->one(FLERR,"Invalid atom ID {} in Velocities section of data file: {}", tagdata, buf);
if ((m = map(tagdata)) >= 0) avec->data_vel(m,values);
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
process N bonds read into buf from data files
if count is non-nullptr, 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 labelflag, int *ilabel)
{
int m,itype;
tagint atom1,atom2;
char *next;
std::string typestr;
int newton_bond = force->newton_bond;
auto location = "Bonds section of data file";
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
if (!next) error->all(FLERR, "Missing data in {}", location);
*next = '\0';
auto values = Tokenizer(buf).as_vector();
int nwords = values.size();
for (std::size_t ii = 0; ii < values.size(); ++ii) {
if (utils::strmatch(values[ii], "^#")) {
nwords = ii;
break;
}
}
// skip over empty or comment lines
// Bonds line is: number(ignored), bond type, atomID 1, atomID 2
if (nwords > 0) {
if (nwords != 4) error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
typestr = utils::utf8_subst(values[1]);
atom1 = utils::tnumeric(FLERR, values[2], false, lmp);
atom2 = utils::tnumeric(FLERR, values[3], false, lmp);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
}
switch (utils::is_type(typestr)) {
case 0: { // numeric
itype = utils::inumeric(FLERR, typestr, false, lmp) + type_offset;
if ((itype < 1) || (itype > nbondtypes))
error->all(FLERR, "Invalid bond type {} in {}: {}", itype, location, utils::trim(buf));
if (labelflag) itype = ilabel[itype - 1];
break;
}
case 1: { // type label
if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
itype = lmap->find(typestr, Atom::BOND);
if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
default: // invalid
error->one(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
if ((atom1 <= 0) || (atom1 > map_tag_max) ||
(atom2 <= 0) || (atom2 > map_tag_max) || (atom1 == atom2))
error->all(FLERR,"Invalid atom ID in {}: {}", location, utils::trim(buf));
if ((itype <= 0) || (itype > nbondtypes))
error->all(FLERR,"Invalid bond type {} in {}: {}", itype, location, utils::trim(buf));
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]++;
avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset);
}
}
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]++;
avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset);
}
}
}
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
process N angles read into buf from data files
if count is non-nullptr, 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 labelflag, int *ilabel)
{
int m,itype;
tagint atom1,atom2,atom3;
char *next;
std::string typestr;
int newton_bond = force->newton_bond;
auto location = "Angles section of data file";
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
if (!next) error->all(FLERR, "Missing data in {}", location);
*next = '\0';
auto values = Tokenizer(buf).as_vector();
int nwords = values.size();
for (std::size_t ii = 0; ii < values.size(); ++ii) {
if (utils::strmatch(values[ii], "^#")) {
nwords = ii;
break;
}
}
// skip over empty or comment lines
// Angles line is: number(ignored), angle type, atomID 1, atomID 2, atomID 3
if (nwords > 0) {
if (nwords != 5) error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
typestr = utils::utf8_subst(values[1]);
atom1 = utils::tnumeric(FLERR, values[2], false, lmp);
atom2 = utils::tnumeric(FLERR, values[3], false, lmp);
atom3 = utils::tnumeric(FLERR, values[4], false, lmp);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
atom3 += id_offset;
}
switch (utils::is_type(typestr)) {
case 0: { // numeric
itype = utils::inumeric(FLERR, typestr, false, lmp) + type_offset;
if ((itype < 1) || (itype > nangletypes))
error->all(FLERR, "Invalid angle type {} in {}: {}", itype, location, utils::trim(buf));
if (labelflag) itype = ilabel[itype - 1];
break;
}
case 1: { // type label
if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
itype = lmap->find(typestr, Atom::ANGLE);
if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
default: // invalid
error->one(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
if ((atom1 <= 0) || (atom1 > map_tag_max) ||
(atom2 <= 0) || (atom2 > map_tag_max) ||
(atom3 <= 0) || (atom3 > map_tag_max) ||
(atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3))
error->one(FLERR,"Invalid atom ID in {}: {}", location, utils::trim(buf));
if (itype <= 0 || itype > nangletypes)
error->one(FLERR,"Invalid angle type {} in {}: {}", itype, location, utils::trim(buf));
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-nullptr, 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 labelflag, int *ilabel)
{
int m,itype;
tagint atom1,atom2,atom3,atom4;
char *next;
std::string typestr;
int newton_bond = force->newton_bond;
auto location = "Dihedrals section of data file";
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
if (!next) error->all(FLERR, "Missing data in {}", location);
*next = '\0';
auto values = Tokenizer(buf).as_vector();
int nwords = values.size();
for (std::size_t ii = 0; ii < values.size(); ++ii) {
if (utils::strmatch(values[ii], "^#")) {
nwords = ii;
break;
}
}
// skip over empty or comment lines
// Dihedrals line is: number(ignored), bond type, atomID 1, atomID 2, atomID 3, atomID 4
if (nwords > 0) {
if (nwords != 6) error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
typestr = utils::utf8_subst(values[1]);
atom1 = utils::tnumeric(FLERR, values[2], false, lmp);
atom2 = utils::tnumeric(FLERR, values[3], false, lmp);
atom3 = utils::tnumeric(FLERR, values[4], false, lmp);
atom4 = utils::tnumeric(FLERR, values[5], false, lmp);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
atom3 += id_offset;
atom4 += id_offset;
}
switch (utils::is_type(typestr)) {
case 0: { // numeric
itype = utils::inumeric(FLERR, typestr, false, lmp) + type_offset;
if ((itype < 1) || (itype > ndihedraltypes))
error->all(FLERR, "Invalid dihedral type {} in {}: {}", itype, location,
utils::trim(buf));
if (labelflag) itype = ilabel[itype - 1];
break;
}
case 1: { // type label
if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
itype = lmap->find(typestr, Atom::DIHEDRAL);
if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
default: // invalid
error->one(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
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) ||
(atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
(atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
error->one(FLERR, "Invalid atom ID in {}: {}", location, utils::trim(buf));
if (itype <= 0 || itype > ndihedraltypes)
error->one(FLERR, "Invalid dihedral type {} in {}: {}", itype, location, utils::trim(buf));
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-nullptr, 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 labelflag, int *ilabel)
{
int m,itype;
tagint atom1,atom2,atom3,atom4;
char *next;
std::string typestr;
int newton_bond = force->newton_bond;
auto location = "Impropers section of data file";
for (int i = 0; i < n; i++) {
next = strchr(buf,'\n');
if (!next) error->all(FLERR, "Missing data in {}", location);
*next = '\0';
auto values = Tokenizer(buf).as_vector();
int nwords = values.size();
for (std::size_t ii = 0; ii < values.size(); ++ii) {
if (utils::strmatch(values[ii], "^#")) {
nwords = ii;
break;
}
}
// skip over empty or comment lines
// Impropers line is: number(ignored), bond type, atomID 1, atomID 2, atomID 3, atomID 4
if (nwords > 0) {
if (nwords != 6) error->all(FLERR, "Incorrect format in {}: {}", location, utils::trim(buf));
typestr = utils::utf8_subst(values[1]);
atom1 = utils::tnumeric(FLERR, values[2], false, lmp);
atom2 = utils::tnumeric(FLERR, values[3], false, lmp);
atom3 = utils::tnumeric(FLERR, values[4], false, lmp);
atom4 = utils::tnumeric(FLERR, values[5], false, lmp);
if (id_offset) {
atom1 += id_offset;
atom2 += id_offset;
atom3 += id_offset;
atom4 += id_offset;
}
switch (utils::is_type(typestr)) {
case 0: { // numeric
itype = utils::inumeric(FLERR, typestr, false, lmp) + type_offset;
if ((itype < 1) || (itype > nimpropertypes))
error->all(FLERR, "Invalid improper type {} in {}: {}", itype, location,
utils::trim(buf));
if (labelflag) itype = ilabel[itype - 1];
break;
}
case 1: { // type label
if (!atom->labelmapflag) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
itype = lmap->find(typestr, Atom::IMPROPER);
if (itype == -1) error->all(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
default: // invalid
error->one(FLERR, "Invalid {}: {}", location, utils::trim(buf));
break;
}
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) ||
(atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
(atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
error->one(FLERR, "Invalid atom ID in {}: {}", location, utils::trim(buf));
if (itype <= 0 || itype > nimpropertypes)
error->one(FLERR, "Invalid improper type {} in {}: {}", itype, location, utils::trim(buf));
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 m;
char *next;
// 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');
if (!next) error->all(FLERR, "Missing data in Bonus section of data file");
*next = '\0';
auto values = Tokenizer(utils::trim_comment(buf)).as_vector();
if (values.size() == 0) {
// skip over empty or comment lines
} else if ((int)values.size() != avec_bonus->size_data_bonus) {
error->all(FLERR, "Incorrect bonus data format in data file: {}", utils::trim(buf));
} else {
tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + 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);
}
buf = next + 1;
}
}
/* ----------------------------------------------------------------------
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, AtomVec *avec_body, tagint id_offset)
{
std::vector<int> ivalues;
std::vector<double> dvalues;
if (!unique_tags) unique_tags = new std::set<tagint>;
// 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++) {
char *next = strchr(buf,'\n');
if (!next) error->all(FLERR, "Missing data in Bodies section of data file");
*next = '\0';
auto values = Tokenizer(utils::trim_comment(buf)).as_vector();
if (values.size()) {
tagint tagdata = utils::tnumeric(FLERR,values[0],false,lmp) + id_offset;
int ninteger = utils::inumeric(FLERR,values[1],false,lmp);
int ndouble = utils::inumeric(FLERR,values[2],false,lmp);
if (unique_tags->find(tagdata) == unique_tags->end())
unique_tags->insert(tagdata);
else
error->one(FLERR,"Duplicate atom ID {} in Bodies section of data file", tagdata);
buf = next + 1;
int m = map(tagdata);
if (m >= 0) {
ivalues.resize(ninteger);
dvalues.resize(ndouble);
for (int j = 0; j < ninteger; j++) {
buf += strspn(buf," \t\n\r\f");
buf[strcspn(buf," \t\n\r\f")] = '\0';
ivalues[j] = utils::inumeric(FLERR,buf,false,lmp);
buf += strlen(buf)+1;
}
for (int j = 0; j < ndouble; j++) {
buf += strspn(buf," \t\n\r\f");
buf[strcspn(buf," \t\n\r\f")] = '\0';
dvalues[j] = utils::numeric(FLERR,buf,false,lmp);
buf += strlen(buf)+1;
}
avec_body->data_body(m,ninteger,ndouble,ivalues.data(),dvalues.data());
} else {
int nvalues = ninteger + ndouble; // number of values to skip
for (int j = 0; j < nvalues; j++) {
buf += strspn(buf," \t\n\r\f");
buf[strcspn(buf," \t\n\r\f")] = '\0';
buf += strlen(buf)+1;
}
}
}
buf += strspn(buf," \t\n\r\f");
}
}
/* ----------------------------------------------------------------------
init per-atom fix/compute/variable values for newly created atoms
called from create_atoms, read_data, read_dump,
lib::lammps_create_atoms()
fixes, computes, variables may or may not exist when called
------------------------------------------------------------------------- */
void Atom::data_fix_compute_variable(int nprev, int nnew)
{
for (const auto &fix : modify->get_fix_list()) {
if (fix->create_attribute)
for (int i = nprev; i < nnew; i++)
fix->set_arrays(i);
}
for (int m = 0; m < modify->ncompute; m++) {
Compute *compute = modify->compute[m];
if (compute->create_attribute)
for (int i = nprev; i < nnew; i++)
compute->set_arrays(i);
}
for (int i = nprev; i < nnew; i++)
input->variable->set_arrays(i);
}
/* ----------------------------------------------------------------------
allocate arrays of length ntypes
only done after ntypes is set
------------------------------------------------------------------------- */
void Atom::allocate_type_arrays()
{
if (avec->mass_type == AtomVec::PER_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
------------------------------------------------------------------------- */
// clang-format on
void Atom::set_mass(const char *file, int line, const char *str, int type_offset, int labelflag,
int *ilabel)
{
if (mass == nullptr) error->all(file, line, "Cannot set mass for atom style {}", atom_style);
int itype;
double mass_one;
auto location = "Masses section of data file";
auto values = Tokenizer(str).as_vector();
int nwords = values.size();
for (std::size_t i = 0; i < values.size(); ++i) {
if (utils::strmatch(values[i], "^#")) {
nwords = i;
break;
}
}
if (nwords != 2) error->all(file, line, "Invalid format in {}: {}", location, str);
auto typestr = utils::utf8_subst(values[0]);
switch (utils::is_type(typestr)) {
case 0: { // numeric
itype = utils::inumeric(file, line, typestr, false, lmp);
if ((itype < 1) || (itype > ntypes))
error->all(file, line, "Invalid atom type {} in {}: {}", itype, location, utils::trim(str));
if (labelflag) itype = ilabel[itype - 1];
break;
}
case 1: { // type label
if (!atom->labelmapflag)
error->all(file, line, "Invalid atom type in {}: {}", location, utils::trim(str));
itype = lmap->find(typestr, Atom::ATOM);
if (itype == -1)
error->all(file, line, "Unknown atom type {} in {}: {}", typestr, location,
utils::trim(str));
break;
}
default: // invalid
itype = -1000000000;
error->one(file, line, "Invalid {}: {}", location, utils::trim(str));
break;
}
itype += type_offset;
mass_one = utils::numeric(file, line, values[1], false, lmp);
if (itype < 1 || itype > ntypes)
error->all(file, line, "Invalid atom type {} in {}: {}", itype, location, utils::trim(str));
if (mass_one <= 0.0)
error->all(file, line, "Invalid mass value {} in {}: {}", mass_one, location, utils::trim(str));
mass[itype] = mass_one;
mass_setflag[itype] = 1;
}
// clang-format off
/* ----------------------------------------------------------------------
set a mass and flag it as set
called from EAM pair routine
------------------------------------------------------------------------- */
void Atom::set_mass(const char *file, int line, int itype, double value)
{
if (mass == nullptr) error->all(file,line, "Cannot set mass for atom style {}", atom_style);
if (itype < 1 || itype > ntypes)
error->all(file,line,"Invalid type {} for atom mass {}", itype, value);
if (value <= 0.0) error->all(file,line,"Invalid atom mass value {}", value);
mass[itype] = value;
mass_setflag[itype] = 1;
}
/* ----------------------------------------------------------------------
set one or more masses and flag them as set
called from reading of input script
------------------------------------------------------------------------- */
void Atom::set_mass(const char *file, int line, int /*narg*/, char **arg)
{
if (mass == nullptr)
error->all(file,line, "Cannot set per-type atom mass for atom style {}", atom_style);
char *typestr = utils::expand_type(file, line, arg[0], Atom::ATOM, lmp);
if (typestr) arg[0] = typestr;
int lo, hi;
utils::bounds(file, line, arg[0], 1, ntypes, lo, hi, error);
if ((lo < 1) || (hi > ntypes))
error->all(file, line, "Invalid atom type {} for atom mass", arg[0]);
const double value = utils::numeric(FLERR, arg[1], false, lmp);
if (value <= 0.0)
error->all(file, line, "Invalid atom mass value {} for type {}", value, arg[0]);
for (int itype = lo; itype <= hi; itype++) {
mass[itype] = value;
mass_setflag[itype] = 1;
}
}
/* ----------------------------------------------------------------------
set all masses
called from reading of restart file, also from ServerMD
------------------------------------------------------------------------- */
void Atom::set_mass(double *values)
{
for (int itype = 1; itype <= ntypes; itype++) {
mass[itype] = values[itype];
mass_setflag[itype] = 1;
}
}
/* ----------------------------------------------------------------------
check that all per-atom-type masses have been set
------------------------------------------------------------------------- */
void Atom::check_mass(const char *file, int line)
{
if (mass == nullptr) return;
if (rmass_flag) return;
for (int itype = 1; itype <= ntypes; itype++)
if (mass_setflag[itype] == 0)
error->all(file,line,"Not all per-type masses are set. Type {} is missing.", itype);
}
/* ----------------------------------------------------------------------
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;
auto avec_ellipsoid = dynamic_cast<AtomVecEllipsoid *>(style_match("ellipsoid"));
auto 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) utils::missing_cmd_args(FLERR, "molecule", error);
if (find_molecule(arg[0]) >= 0)
error->all(FLERR,"Reuse of molecule template ID {}", arg[0]);
// 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 (true) {
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(const char *id)
{
if (id == nullptr) return -1;
for (int imol = 0; imol < nmolecule; imol++)
if (strcmp(id,molecules[imol]->id) == 0) return imol;
return -1;
}
/* ----------------------------------------------------------------------
return vector of molecules which match template ID
------------------------------------------------------------------------- */
std::vector<Molecule *>Atom::get_molecule_by_id(const std::string &id)
{
std::vector<Molecule *> result;
for (int imol = 0; imol < nmolecule; ++imol)
if (id == molecules[imol]->id) result.push_back(molecules[imol]);
return result;
}
/* ----------------------------------------------------------------------
add info to current atom ilocal from molecule template onemol and its iatom
offset = atom ID preceding 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);
}
// initialize custom per-atom properties to zero if present
for (int i = 0; i < nivector; ++i) ivector[i][ilocal] = 0;
for (int i = 0; i < ndvector; ++i) dvector[i][ilocal] = 0.0;
for (int i = 0; i < niarray; ++i)
for (int j = 0; j < icols[i]; ++j)
iarray[i][ilocal][j] = 0;
for (int i = 0; i < ndarray; ++i)
for (int j = 0; j < dcols[i]; ++j)
darray[i][ilocal][j] = 0.0;
if (molecular != Atom::MOLECULAR) 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;
}
}
/* ----------------------------------------------------------------------
allocate space for type label map
------------------------------------------------------------------------- */
void Atom::add_label_map()
{
if (lmp->kokkos)
error->all(FLERR, "Label maps are currently not supported with Kokkos");
labelmapflag = 1;
lmap = new LabelMap(lmp,ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes);
}
/* ----------------------------------------------------------------------
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;
// 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];
}
// 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
// default = 1/2 of neighbor cutoff
// check if neighbor cutoff = 0.0
// and in that case, disable sorting
double binsize = 0.0;
if (userbinsize > 0.0) binsize = userbinsize;
else if (neighbor->cutneighmax > 0.0) binsize = 0.5 * neighbor->cutneighmax;
if ((binsize == 0.0) && (sortfreq > 0)) {
sortfreq = 0;
if (comm->me == 0)
error->warning(FLERR,"No pairwise cutoff or binsize set. Atom sorting therefore disabled.");
return;
}
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]);
#ifdef LMP_INTEL
if (neighbor->has_intel_request() && userbinsize == 0.0) {
if (neighbor->binsizeflag) bininv = 1.0/neighbor->binsize_user;
double nx_low = neighbor->bboxlo[0];
double ny_low = neighbor->bboxlo[1];
double nz_low = neighbor->bboxlo[2];
double nxbbox = neighbor->bboxhi[0] - nx_low;
double nybbox = neighbor->bboxhi[1] - ny_low;
double nzbbox = neighbor->bboxhi[2] - nz_low;
int nnbinx = static_cast<int> (nxbbox * bininv);
int nnbiny = static_cast<int> (nybbox * bininv);
int nnbinz = static_cast<int> (nzbbox * bininv);
if (domain->dimension == 2) nnbinz = 1;
if (nnbinx == 0) nnbinx = 1;
if (nnbiny == 0) nnbiny = 1;
if (nnbinz == 0) nnbinz = 1;
double binsizex = nxbbox/nnbinx;
double binsizey = nybbox/nnbiny;
double binsizez = nzbbox/nnbinz;
bininvx = 1.0 / binsizex;
bininvy = 1.0 / binsizey;
bininvz = 1.0 / binsizez;
int lxo = (bboxlo[0] - nx_low) * bininvx;
int lyo = (bboxlo[1] - ny_low) * bininvy;
int lzo = (bboxlo[2] - nz_low) * bininvz;
bboxlo[0] = nx_low + static_cast<double>(lxo) / bininvx;
bboxlo[1] = ny_low + static_cast<double>(lyo) / bininvy;
bboxlo[2] = nz_low + static_cast<double>(lzo) / bininvz;
nbinx = static_cast<int>((bboxhi[0] - bboxlo[0]) * bininvx) + 1;
nbiny = static_cast<int>((bboxhi[1] - bboxlo[1]) * bininvy) + 1;
nbinz = static_cast<int>((bboxhi[2] - bboxlo[2]) * bininvz) + 1;
bboxhi[0] = bboxlo[0] + static_cast<double>(nbinx) / bininvx;
bboxhi[1] = bboxlo[1] + static_cast<double>(nbiny) / bininvy;
bboxhi[2] = bboxlo[2] + static_cast<double>(nbinz) / bininvz;
}
#endif
#ifdef LMP_GPU
if (userbinsize == 0.0) {
FixGPU *fix = dynamic_cast<FixGPU *>(modify->get_fix_by_id("package_gpu"));
if (fix) {
const double subx = domain->subhi[0] - domain->sublo[0];
const double suby = domain->subhi[1] - domain->sublo[1];
const double subz = domain->subhi[2] - domain->sublo[2];
binsize = fix->binsize(subx, suby, subz, atom->nlocal,neighbor->cutneighmax);
bininv = 1.0 / binsize;
nbinx = static_cast<int> (ceil(subx * bininv));
nbiny = static_cast<int> (ceil(suby * bininv));
nbinz = static_cast<int> (ceil(subz * bininv));
if (domain->dimension == 2) nbinz = 1;
if (nbinx == 0) nbinx = 1;
if (nbiny == 0) nbiny = 1;
if (nbinz == 0) nbinz = 1;
bininvx = bininv;
bininvy = bininv;
bininvz = bininv;
}
}
#endif
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 = Atom::GROW for grow, Atom::RESTART for restart, Atom::BORDER for border comm
------------------------------------------------------------------------- */
void Atom::add_callback(int flag)
{
int ifix;
// find the fix
// if find null pointer:
// 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 pointer:
// 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] == nullptr) break;
// add callback to lists and sort, reallocating if necessary
// sorting is required in cases where fixes were replaced as it ensures atom
// data is read/written/transfered in the same order that fixes are called
if (flag == GROW) {
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++;
std::sort(extra_grow, extra_grow + nextra_grow);
} else if (flag == RESTART) {
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++;
std::sort(extra_restart, extra_restart + nextra_restart);
} else if (flag == BORDER) {
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++;
std::sort(extra_border, extra_border + 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)
{
if (id == nullptr) return;
int ifix = modify->find_fix(id);
// compact the list of callbacks
if (flag == GROW) {
int match;
for (match = 0; match < nextra_grow; match++)
if (extra_grow[match] == ifix) break;
if ((nextra_grow == 0) || (match == nextra_grow))
error->all(FLERR,"Trying to delete non-existent Atom::grow() callback");
for (int i = match; i < nextra_grow-1; i++)
extra_grow[i] = extra_grow[i+1];
nextra_grow--;
} else if (flag == RESTART) {
int match;
for (match = 0; match < nextra_restart; match++)
if (extra_restart[match] == ifix) break;
if ((nextra_restart == 0) || (match == nextra_restart))
error->all(FLERR,"Trying to delete non-existent Atom::restart() callback");
for (int i = match; i < nextra_restart-1; i++)
extra_restart[i] = extra_restart[i+1];
nextra_restart--;
} else if (flag == BORDER) {
int match;
for (match = 0; match < nextra_border; match++)
if (extra_border[match] == ifix) break;
if ((nextra_border == 0) || (match == nextra_border))
error->all(FLERR,"Trying to delete non-existent Atom::border() callback");
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, -1 if not found
lists of names can have NULL entries if previously removed
return flag = 0/1 for int/double
return cols = 0/N for vector/array where N = # of columns
------------------------------------------------------------------------- */
int Atom::find_custom(const char *name, int &flag, int &cols)
{
if (name == nullptr) return -1;
for (int i = 0; i < nivector; i++)
if (ivname[i] && strcmp(ivname[i],name) == 0) {
flag = 0;
cols = 0;
return i;
}
for (int i = 0; i < ndvector; i++)
if (dvname[i] && strcmp(dvname[i],name) == 0) {
flag = 1;
cols = 0;
return i;
}
for (int i = 0; i < niarray; i++)
if (ianame[i] && strcmp(ianame[i],name) == 0) {
flag = 0;
cols = icols[i];
return i;
}
for (int i = 0; i < ndarray; i++)
if (daname[i] && strcmp(daname[i],name) == 0) {
flag = 1;
cols = dcols[i];
return i;
}
return -1;
}
/** \brief Add a custom per-atom property with the given name and type and size
\verbatim embed:rst
This function will add a custom per-atom property with one or more values
with the name "name" to the list of custom properties.
This function is called, e.g. from :doc:`fix property/atom <fix_property_atom>`.
\endverbatim
* \param name Name of the property (w/o a "i_" or "d_" or "i2_" or "d2_" prefix)
* \param flag Data type of property: 0 for int, 1 for double
* \param cols Number of values: 0 for a single value, 1 or more for a vector of values
* \return index of property in the respective list of properties
*/
int Atom::add_custom(const char *name, int flag, int cols)
{
int index = -1;
if ((flag == 0) && (cols == 0)) {
index = nivector;
nivector++;
ivname = (char **) memory->srealloc(ivname,nivector*sizeof(char *),"atom:ivname");
ivname[index] = utils::strdup(name);
ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),"atom:ivector");
memory->create(ivector[index],nmax,"atom:ivector");
} else if ((flag == 1) && (cols == 0)) {
index = ndvector;
ndvector++;
dvname = (char **) memory->srealloc(dvname,ndvector*sizeof(char *),"atom:dvname");
dvname[index] = utils::strdup(name);
dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),"atom:dvector");
memory->create(dvector[index],nmax,"atom:dvector");
} else if ((flag == 0) && (cols > 0)) {
index = niarray;
niarray++;
ianame = (char **) memory->srealloc(ianame,niarray*sizeof(char *),"atom:ianame");
ianame[index] = utils::strdup(name);
iarray = (int ***) memory->srealloc(iarray,niarray*sizeof(int **),"atom:iarray");
memory->create(iarray[index],nmax,cols,"atom:iarray");
icols = (int *) memory->srealloc(icols,niarray*sizeof(int),"atom:icols");
icols[index] = cols;
} else if ((flag == 1) && (cols > 0)) {
index = ndarray;
ndarray++;
daname = (char **) memory->srealloc(daname,ndarray*sizeof(char *),"atom:daname");
daname[index] = utils::strdup(name);
darray = (double ***) memory->srealloc(darray,ndarray*sizeof(double **),"atom:darray");
memory->create(darray[index],nmax,cols,"atom:darray");
dcols = (int *) memory->srealloc(dcols,ndarray*sizeof(int),"atom:dcols");
dcols[index] = cols;
}
if (index < 0)
error->all(FLERR,"Invalid call to Atom::add_custom()");
return index;
}
/*! \brief Remove a custom per-atom property of a given type and size
*
\verbatim embed:rst
This will remove a property that was requested, e.g. by the
:doc:`fix property/atom <fix_property_atom>` command. It frees the
allocated memory and sets the pointer to ``nullptr`` for the entry in
the list so it can be reused. The lists of these pointers are never
compacted or shrink, so that indices to name mappings remain valid.
\endverbatim
* \param index Index of property in the respective list of properties
* \param flag Data type of property: 0 for int, 1 for double
* \param cols Number of values: 0 for a single value, 1 or more for a vector of values
*/
void Atom::remove_custom(int index, int flag, int cols)
{
if (flag == 0 && cols == 0) {
memory->destroy(ivector[index]);
ivector[index] = nullptr;
delete[] ivname[index];
ivname[index] = nullptr;
} else if (flag == 1 && cols == 0) {
memory->destroy(dvector[index]);
dvector[index] = nullptr;
delete[] dvname[index];
dvname[index] = nullptr;
} else if (flag == 0 && cols) {
memory->destroy(iarray[index]);
iarray[index] = nullptr;
delete[] ianame[index];
ianame[index] = nullptr;
} else if (flag == 1 && cols) {
memory->destroy(darray[index]);
darray[index] = nullptr;
delete[] daname[index];
daname[index] = nullptr;
}
}
/** Provide access to internal data of the Atom class by keyword
*
\verbatim embed:rst
This function is a way to access internal per-atom data. This data is
distributed across MPI ranks and thus only the data for "local" atoms
are expected to be available. Whether also data for "ghost" atoms is
stored and up-to-date depends on various simulation settings.
This table lists a large part of the supported names, their data types,
length of the data area, and a short description.
.. list-table::
:header-rows: 1
:widths: auto
* - Name
- Type
- Items per atom
- Description
* - mass
- double
- 1
- per-type mass. This array is **NOT** a per-atom array
but of length ``ntypes+1``, element 0 is ignored.
* - id
- tagint
- 1
- atom ID of the particles
* - type
- int
- 1
- atom type of the particles
* - mask
- int
- 1
- bitmask for mapping to groups. Individual bits are set
to 0 or 1 for each group.
* - image
- imageint
- 1
- 3 image flags encoded into a single integer.
See :cpp:func:`lammps_encode_image_flags`.
* - x
- double
- 3
- x-, y-, and z-coordinate of the particles
* - v
- double
- 3
- x-, y-, and z-component of the velocity of the particles
* - f
- double
- 3
- x-, y-, and z-component of the force on the particles
* - molecule
- int
- 1
- molecule ID of the particles
* - q
- double
- 1
- charge of the particles
* - mu
- double
- 3
- dipole moment of the particles
* - omega
- double
- 3
- x-, y-, and z-component of rotational velocity of the particles
* - angmom
- double
- 3
- x-, y-, and z-component of angular momentum of the particles
* - torque
- double
- 3
- x-, y-, and z-component of the torque on the particles
* - radius
- double
- 1
- radius of the (extended) particles
* - rmass
- double
- 1
- per-atom mass of the particles. ``nullptr`` if per-type masses are
used. See the :cpp:func:`rmass_flag<lammps_extract_setting>` setting.
* - ellipsoid
- int
- 1
- 1 if the particle is an ellipsoidal particle, 0 if not
* - line
- int
- 1
- 1 if the particle is a line particle, 0 if not
* - tri
- int
- 1
- 1 if the particle is a triangulated particle, 0 if not
* - body
- int
- 1
- 1 if the particle is a body particle, 0 if not
* - quat
- double
- 4
- four quaternion components of the particles
* - i_name
- int
- 1
- single integer value defined by fix property/atom vector name
* - d_name
- double
- 1
- single double value defined by fix property/atom vector name
* - i2_name
- int
- n
- N integer values defined by fix property/atom array name
* - d2_name
- double
- n
- N double values defined by fix property/atom array name
*See also*
:cpp:func:`lammps_extract_atom`
\endverbatim
*
* \sa extract_datatype
*
* \param name string with the keyword of the desired property.
Typically the name of the pointer variable returned
* \return pointer to the requested data cast to ``void *`` or ``nullptr`` */
void *Atom::extract(const char *name)
{
// --------------------------------------------------------------------
// 4th customization section: customize by adding new variable name
// if new variable is from a package, add package comment
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,"body") == 0) return (void *) body;
if (strcmp(name,"quat") == 0) return (void *) quat;
// PERI PACKAGE
if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
if (strcmp(name,"s0") == 0) return (void *) s0;
if (strcmp(name,"x0") == 0) return (void *) x0;
// SPIN PACKAGE
if (strcmp(name,"sp") == 0) return (void *) sp;
// EFF and AWPMD packages
if (strcmp(name,"espin") == 0) return (void *) spin;
if (strcmp(name,"spin") == 0) return (void *) spin; // backward compatibility
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;
// SPH package
if (strcmp(name,"rho") == 0) return (void *) rho;
if (strcmp(name,"drho") == 0) return (void *) drho;
if (strcmp(name,"esph") == 0) return (void *) esph;
if (strcmp(name,"desph") == 0) return (void *) desph;
if (strcmp(name,"cv") == 0) return (void *) cv;
if (strcmp(name,"vest") == 0) return (void *) vest;
// MESONT package
if (strcmp(name,"length") == 0) return (void *) length;
if (strcmp(name,"buckling") == 0) return (void *) buckling;
if (strcmp(name,"bond_nt") == 0) return (void *) bond_nt;
// MACHDYN package
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;
// DPD-REACT pakage
if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
// DPD-MESO package
if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp;
// DIELECTRIC package
if (strcmp(name,"area") == 0) return (void *) area;
if (strcmp(name,"ed") == 0) return (void *) ed;
if (strcmp(name,"em") == 0) return (void *) em;
if (strcmp(name,"epsilon") == 0) return (void *) epsilon;
if (strcmp(name,"curvature") == 0) return (void *) curvature;
if (strcmp(name,"q_unscaled") == 0) return (void *) q_unscaled;
// end of customization section
// --------------------------------------------------------------------
// custom vectors and arrays
if (utils::strmatch(name,"^[id]2?_")) {
int which = 0, array = 0;
if (name[0] == 'd') which = 1;
if (name[1] == '2') array = 1;
int index,flag,cols;
if (!array) index = find_custom(&name[2],flag,cols);
else index = find_custom(&name[3],flag,cols);
if (index < 0) return nullptr;
if (which != flag) return nullptr;
if ((!array && cols) || (array && !cols)) return nullptr;
if (!which && !array) return (void *) ivector[index];
if (which && !array) return (void *) dvector[index];
if (!which && array) return (void *) iarray[index];
if (which && array) return (void *) darray[index];
}
return nullptr;
}
/** Provide data type info about internal data of the Atom class
*
\verbatim embed:rst
.. versionadded:: 18Sep2020
\endverbatim
*
* \sa extract
*
* \param name string with the keyword of the desired property.
* \return data type constant for desired property or -1 */
int Atom::extract_datatype(const char *name)
{
// --------------------------------------------------------------------
// 5th customization section: customize by adding new variable name
if (strcmp(name,"mass") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"id") == 0) return LAMMPS_TAGINT;
if (strcmp(name,"type") == 0) return LAMMPS_INT;
if (strcmp(name,"mask") == 0) return LAMMPS_INT;
if (strcmp(name,"image") == 0) return LAMMPS_TAGINT;
if (strcmp(name,"x") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"v") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"f") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"molecule") == 0) return LAMMPS_TAGINT;
if (strcmp(name,"q") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"mu") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"omega") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"angmom") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"torque") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"radius") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"rmass") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"ellipsoid") == 0) return LAMMPS_INT;
if (strcmp(name,"line") == 0) return LAMMPS_INT;
if (strcmp(name,"tri") == 0) return LAMMPS_INT;
if (strcmp(name,"body") == 0) return LAMMPS_INT;
if (strcmp(name,"quat") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"x0") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"espin") == 0) return LAMMPS_INT;
if (strcmp(name,"spin") == 0) return LAMMPS_INT; // backwards compatibility
if (strcmp(name,"eradius") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"ervel") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"erforce") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"ervelforce") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"cs") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"csforce") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"vforce") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name,"etag") == 0) return LAMMPS_INT;
if (strcmp(name,"rho") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"drho") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"esph") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"desph") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"cv") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"vest") == 0) return LAMMPS_DOUBLE_2D;
// MESONT package
if (strcmp(name,"length") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"buckling") == 0) return LAMMPS_INT;
if (strcmp(name,"bond_nt") == 0) return LAMMPS_TAGINT_2D;
// MACHDYN package
if (strcmp(name, "contact_radius") == 0) return LAMMPS_DOUBLE;
if (strcmp(name, "smd_data_9") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name, "smd_stress") == 0) return LAMMPS_DOUBLE_2D;
if (strcmp(name, "eff_plastic_strain") == 0) return LAMMPS_DOUBLE;
if (strcmp(name, "eff_plastic_strain_rate") == 0) return LAMMPS_DOUBLE;
if (strcmp(name, "damage") == 0) return LAMMPS_DOUBLE;
// DPD-REACT package
if (strcmp(name,"dpdTheta") == 0) return LAMMPS_DOUBLE;
// DPD-MESO package
if (strcmp(name,"edpd_temp") == 0) return LAMMPS_DOUBLE;
// DIELECTRIC package
if (strcmp(name,"area") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"ed") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"em") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"epsilon") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"curvature") == 0) return LAMMPS_DOUBLE;
if (strcmp(name,"q_unscaled") == 0) return LAMMPS_DOUBLE;
// end of customization section
// --------------------------------------------------------------------
// custom vectors and arrays
if (utils::strmatch(name,"^[id]2?_")) {
int which = 0, array = 0;
if (name[0] == 'd') which = 1;
if (name[1] == '2') array = 1;
int index,flag,cols;
if (!array) index = find_custom(&name[2],flag,cols);
else index = find_custom(&name[3],flag,cols);
if (index < 0) return -1;
if (which != flag) return -1;
if ((!array && cols) || (array && !cols)) return -1;
if (which == 0) return LAMMPS_INT;
else return LAMMPS_DOUBLE;
}
return -1;
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
call to avec tallies per-atom vectors
add in global to local mapping storage
------------------------------------------------------------------------- */
double Atom::memory_usage()
{
double bytes = avec->memory_usage();
bytes += (double)max_same*sizeof(int);
if (map_style == MAP_ARRAY)
bytes += memory->usage(map_array,map_maxarray);
else if (map_style == MAP_HASH) {
bytes += (double)map_nbucket*sizeof(int);
bytes += (double)map_nhash*sizeof(HashElem);
}
if (maxnext) {
bytes += memory->usage(next,maxnext);
bytes += memory->usage(permute,maxnext);
}
return bytes;
}
Reference in New Issue View Git Blame Copy Permalink