forked from lijiext/lammps
685 lines
22 KiB
C++
685 lines
22 KiB
C++
/* ----------------------------------------------------------------------
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
|
http://lammps.sandia.gov, Sandia National Laboratories
|
|
Steve Plimpton, sjplimp@sandia.gov
|
|
|
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
|
certain rights in this software. This software is distributed under
|
|
the GNU General Public License.
|
|
|
|
See the README file in the top-level LAMMPS directory.
|
|
------------------------------------------------------------------------- */
|
|
|
|
#include "mpi.h"
|
|
#include "math.h"
|
|
#include "stdlib.h"
|
|
#include "string.h"
|
|
#include "set.h"
|
|
#include "atom.h"
|
|
#include "atom_vec.h"
|
|
#include "atom_vec_ellipsoid.h"
|
|
#include "domain.h"
|
|
#include "region.h"
|
|
#include "group.h"
|
|
#include "comm.h"
|
|
#include "neighbor.h"
|
|
#include "force.h"
|
|
#include "pair.h"
|
|
#include "random_park.h"
|
|
#include "math_extra.h"
|
|
#include "error.h"
|
|
|
|
using namespace LAMMPS_NS;
|
|
|
|
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
|
|
enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,
|
|
DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,
|
|
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
|
|
MESO_E,MESO_CV,MESO_RHO};
|
|
|
|
#define BIG INT_MAX
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
Set::Set(LAMMPS *lmp) : Pointers(lmp)
|
|
{
|
|
PI = 4.0*atan(1.0);
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void Set::command(int narg, char **arg)
|
|
{
|
|
if (domain->box_exist == 0)
|
|
error->all(FLERR,"Set command before simulation box is defined");
|
|
if (atom->natoms == 0)
|
|
error->all(FLERR,"Set command with no atoms existing");
|
|
if (narg < 3) error->all(FLERR,"Illegal set command");
|
|
|
|
// style and ID info
|
|
|
|
if (strcmp(arg[0],"atom") == 0) style = ATOM_SELECT;
|
|
else if (strcmp(arg[0],"mol") == 0) style = MOL_SELECT;
|
|
else if (strcmp(arg[0],"type") == 0) style = TYPE_SELECT;
|
|
else if (strcmp(arg[0],"group") == 0) style = GROUP_SELECT;
|
|
else if (strcmp(arg[0],"region") == 0) style = REGION_SELECT;
|
|
else error->all(FLERR,"Illegal set command");
|
|
|
|
int n = strlen(arg[1]) + 1;
|
|
id = new char[n];
|
|
strcpy(id,arg[1]);
|
|
select = NULL;
|
|
|
|
// loop over keyword/value pairs
|
|
// call appropriate routine to reset attributes
|
|
|
|
if (comm->me == 0 && screen) fprintf(screen,"Setting atom values ...\n");
|
|
|
|
int allcount,origarg;
|
|
|
|
int iarg = 2;
|
|
while (iarg < narg) {
|
|
count = 0;
|
|
origarg = iarg;
|
|
|
|
if (strcmp(arg[iarg],"type") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (ivalue <= 0 || ivalue > atom->ntypes)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
set(TYPE);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"type/fraction") == 0) {
|
|
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
|
|
newtype = atoi(arg[iarg+1]);
|
|
fraction = atof(arg[iarg+2]);
|
|
ivalue = atoi(arg[iarg+3]);
|
|
if (newtype <= 0 || newtype > atom->ntypes)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
if (fraction < 0.0 || fraction > 1.0)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
if (ivalue <= 0) error->all(FLERR,"Invalid random number seed in set command");
|
|
setrandom(TYPE_FRACTION);
|
|
iarg += 4;
|
|
} else if (strcmp(arg[iarg],"mol") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (!atom->molecule_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(MOLECULE);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"x") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
set(X);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"y") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
set(Y);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"z") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
set(Z);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"charge") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->q_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(CHARGE);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"mass") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->rmass_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (dvalue <= 0.0) error->all(FLERR,"Invalid mass in set command");
|
|
set(MASS);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"shape") == 0) {
|
|
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
|
|
xvalue = atof(arg[iarg+1]);
|
|
yvalue = atof(arg[iarg+2]);
|
|
zvalue = atof(arg[iarg+3]);
|
|
if (!atom->ellipsoid_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (xvalue < 0.0 || yvalue < 0.0 || zvalue < 0.0)
|
|
error->all(FLERR,"Invalid shape in set command");
|
|
if (xvalue > 0.0 || yvalue > 0.0 || zvalue > 0.0) {
|
|
if (xvalue == 0.0 || yvalue == 0.0 || zvalue == 0.0)
|
|
error->one(FLERR,"Invalid shape in set command");
|
|
}
|
|
set(SHAPE);
|
|
iarg += 4;
|
|
} else if (strcmp(arg[iarg],"dipole") == 0) {
|
|
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
|
|
xvalue = atof(arg[iarg+1]);
|
|
yvalue = atof(arg[iarg+2]);
|
|
zvalue = atof(arg[iarg+3]);
|
|
if (!atom->mu_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(DIPOLE);
|
|
iarg += 4;
|
|
} else if (strcmp(arg[iarg],"dipole/random") == 0) {
|
|
if (iarg+3 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
dvalue = atof(arg[iarg+2]);
|
|
if (!atom->mu_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (ivalue <= 0) error->all(FLERR,"Invalid random number seed in set command");
|
|
if (dvalue <= 0.0) error->all(FLERR,"Invalid dipole length in set command");
|
|
setrandom(DIPOLE_RANDOM);
|
|
iarg += 3;
|
|
} else if (strcmp(arg[iarg],"quat") == 0) {
|
|
if (iarg+5 > narg) error->all(FLERR,"Illegal set command");
|
|
xvalue = atof(arg[iarg+1]);
|
|
yvalue = atof(arg[iarg+2]);
|
|
zvalue = atof(arg[iarg+3]);
|
|
wvalue = atof(arg[iarg+4]);
|
|
if (!atom->ellipsoid_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(QUAT);
|
|
iarg += 5;
|
|
} else if (strcmp(arg[iarg],"quat/random") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (!atom->ellipsoid_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (ivalue <= 0) error->all(FLERR,"Invalid random number seed in set command");
|
|
setrandom(QUAT_RANDOM);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"diameter") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->radius_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (dvalue < 0.0) error->all(FLERR,"Invalid diameter in set command");
|
|
set(DIAMETER);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"density") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->rmass_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(DENSITY);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"volume") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->vfrac_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(VOLUME);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"image") == 0) {
|
|
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
|
|
ximageflag = yimageflag = zimageflag = 0;
|
|
if (strcmp(arg[iarg+1],"NULL") != 0) {
|
|
ximageflag = 1;
|
|
ximage = atoi(arg[iarg+1]);
|
|
}
|
|
if (strcmp(arg[iarg+2],"NULL") != 0) {
|
|
yimageflag = 1;
|
|
yimage = atoi(arg[iarg+2]);
|
|
}
|
|
if (strcmp(arg[iarg+3],"NULL") != 0) {
|
|
zimageflag = 1;
|
|
zimage = atoi(arg[iarg+3]);
|
|
}
|
|
if (ximageflag && ximage && !domain->xperiodic)
|
|
error->all(FLERR,"Cannot set non-zero image flag for non-periodic dimension");
|
|
if (yimageflag && yimage && !domain->yperiodic)
|
|
error->all(FLERR,"Cannot set non-zero image flag for non-periodic dimension");
|
|
if (zimageflag && zimage && !domain->zperiodic)
|
|
error->all(FLERR,"Cannot set non-zero image flag for non-periodic dimension");
|
|
set(IMAGE);
|
|
iarg += 4;
|
|
} else if (strcmp(arg[iarg],"bond") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (atom->avec->bonds_allow == 0)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (ivalue <= 0 || ivalue > atom->nbondtypes)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
topology(BOND);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"angle") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (atom->avec->angles_allow == 0)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (ivalue <= 0 || ivalue > atom->nangletypes)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
topology(ANGLE);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"dihedral") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (atom->avec->dihedrals_allow == 0)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (ivalue <= 0 || ivalue > atom->ndihedraltypes)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
topology(DIHEDRAL);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"improper") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
ivalue = atoi(arg[iarg+1]);
|
|
if (atom->avec->impropers_allow == 0)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
if (ivalue <= 0 || ivalue > atom->nimpropertypes)
|
|
error->all(FLERR,"Invalid value in set command");
|
|
topology(IMPROPER);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"meso_e") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->e_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(MESO_E);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"meso_cv") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->cv_flag)
|
|
error->all(FLERR,"Cannot set this attribute for this atom style");
|
|
set(MESO_CV);
|
|
iarg += 2;
|
|
} else if (strcmp(arg[iarg],"meso_rho") == 0) {
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
|
|
dvalue = atof(arg[iarg+1]);
|
|
if (!atom->rho_flag)
|
|
error->all(FLERR,"Cannot set meso_rho for this atom style");
|
|
set(MESO_RHO);
|
|
iarg += 2;
|
|
} else error->all(FLERR,"Illegal set command");
|
|
|
|
// statistics
|
|
|
|
MPI_Allreduce(&count,&allcount,1,MPI_INT,MPI_SUM,world);
|
|
|
|
if (comm->me == 0) {
|
|
if (screen) fprintf(screen," %d settings made for %s\n",
|
|
allcount,arg[origarg]);
|
|
if (logfile) fprintf(logfile," %d settings made for %s\n",
|
|
allcount,arg[origarg]);
|
|
}
|
|
}
|
|
|
|
// free local memory
|
|
|
|
delete [] id;
|
|
delete [] select;
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
select atoms according to ATOM, MOLECULE, TYPE, GROUP, REGION style
|
|
n = nlocal or nlocal+nghost depending on keyword
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Set::selection(int n)
|
|
{
|
|
delete [] select;
|
|
select = new int[n];
|
|
int nlo,nhi;
|
|
|
|
if (style == ATOM_SELECT) {
|
|
if (atom->tag_enable == 0)
|
|
error->all(FLERR,"Cannot use set atom with no atom IDs defined");
|
|
force->bounds(id,BIG,nlo,nhi);
|
|
|
|
int *tag = atom->tag;
|
|
for (int i = 0; i < n; i++)
|
|
if (tag[i] >= nlo && tag[i] <= nhi) select[i] = 1;
|
|
else select[i] = 0;
|
|
|
|
} else if (style == MOL_SELECT) {
|
|
if (atom->molecule_flag == 0)
|
|
error->all(FLERR,"Cannot use set mol with no molecule IDs defined");
|
|
if (strcmp(id,"0") == 0) nlo = nhi = 0;
|
|
else force->bounds(id,BIG,nlo,nhi);
|
|
|
|
int *molecule = atom->molecule;
|
|
for (int i = 0; i < n; i++)
|
|
if (molecule[i] >= nlo && molecule[i] <= nhi) select[i] = 1;
|
|
else select[i] = 0;
|
|
|
|
} else if (style == TYPE_SELECT) {
|
|
force->bounds(id,atom->ntypes,nlo,nhi);
|
|
|
|
int *type = atom->type;
|
|
for (int i = 0; i < n; i++)
|
|
if (type[i] >= nlo && type[i] <= nhi) select[i] = 1;
|
|
else select[i] = 0;
|
|
|
|
} else if (style == GROUP_SELECT) {
|
|
int igroup = group->find(id);
|
|
if (igroup == -1) error->all(FLERR,"Could not find set group ID");
|
|
int groupbit = group->bitmask[igroup];
|
|
|
|
int *mask = atom->mask;
|
|
for (int i = 0; i < n; i++)
|
|
if (mask[i] & groupbit) select[i] = 1;
|
|
else select[i] = 0;
|
|
|
|
} else if (style == REGION_SELECT) {
|
|
int iregion = domain->find_region(id);
|
|
if (iregion == -1) error->all(FLERR,"Set region ID does not exist");
|
|
|
|
double **x = atom->x;
|
|
for (int i = 0; i < n; i++)
|
|
if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
|
|
select[i] = 1;
|
|
else select[i] = 0;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
set an owned atom property directly
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Set::set(int keyword)
|
|
{
|
|
AtomVecEllipsoid *avec_ellipsoid =
|
|
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
|
|
|
selection(atom->nlocal);
|
|
|
|
int nlocal = atom->nlocal;
|
|
for (int i = 0; i < nlocal; i++) {
|
|
if (!select[i]) continue;
|
|
|
|
if (keyword == TYPE) atom->type[i] = ivalue;
|
|
else if (keyword == MOLECULE) atom->molecule[i] = ivalue;
|
|
else if (keyword == X) atom->x[i][0] = dvalue;
|
|
else if (keyword == Y) atom->x[i][1] = dvalue;
|
|
else if (keyword == Z) atom->x[i][2] = dvalue;
|
|
else if (keyword == CHARGE) atom->q[i] = dvalue;
|
|
else if (keyword == MASS) atom->rmass[i] = dvalue;
|
|
else if (keyword == DIAMETER) atom->radius[i] = 0.5 * dvalue;
|
|
else if (keyword == VOLUME) atom->vfrac[i] = dvalue;
|
|
else if (keyword == MESO_E) atom->e[i] = dvalue;
|
|
else if (keyword == MESO_CV) atom->cv[i] = dvalue;
|
|
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
|
|
|
|
// set shape
|
|
|
|
else if (keyword == SHAPE)
|
|
avec_ellipsoid->set_shape(i,0.5*xvalue,0.5*yvalue,0.5*zvalue);
|
|
|
|
// set rmass via density
|
|
// if radius > 0.0, treat as sphere
|
|
// if shape > 0.0, treat as ellipsoid
|
|
// else set rmass to density directly
|
|
|
|
else if (keyword == DENSITY) {
|
|
if (atom->radius_flag && atom->radius[i] > 0.0)
|
|
atom->rmass[i] = 4.0*PI/3.0 *
|
|
atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
|
|
else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) {
|
|
double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape;
|
|
atom->rmass[i] = 4.0*PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue;
|
|
} else atom->rmass[i] = dvalue;
|
|
|
|
// reset any or all of 3 image flags
|
|
|
|
} else if (keyword == IMAGE) {
|
|
int xbox = (atom->image[i] & 1023) - 512;
|
|
int ybox = (atom->image[i] >> 10 & 1023) - 512;
|
|
int zbox = (atom->image[i] >> 20) - 512;
|
|
if (ximageflag) xbox = ximage;
|
|
if (yimageflag) ybox = yimage;
|
|
if (zimageflag) zbox = zimage;
|
|
atom->image[i] = ((zbox + 512 & 1023) << 20) |
|
|
((ybox + 512 & 1023) << 10) | (xbox + 512 & 1023);
|
|
|
|
// set dipole moment
|
|
|
|
} else if (keyword == DIPOLE) {
|
|
double **mu = atom->mu;
|
|
mu[i][0] = xvalue;
|
|
mu[i][1] = yvalue;
|
|
mu[i][2] = zvalue;
|
|
mu[i][3] = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] +
|
|
mu[i][2]*mu[i][2]);
|
|
|
|
// set quaternion orientation
|
|
|
|
} else if (keyword == QUAT) {
|
|
if (atom->ellipsoid[i] < 0)
|
|
error->one(FLERR,"Cannot set quaternion for atom that is not an ellipsoid");
|
|
double *quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
|
|
double theta2 = 0.5 * PI * wvalue/180.0;
|
|
double sintheta2 = sin(theta2);
|
|
quat[0] = cos(theta2);
|
|
quat[1] = xvalue * sintheta2;
|
|
quat[2] = yvalue * sintheta2;
|
|
quat[3] = zvalue * sintheta2;
|
|
MathExtra::qnormalize(quat);
|
|
}
|
|
count++;
|
|
}
|
|
}
|
|
|
|
/* ----------------------------------------------------------------------
|
|
set an owned atom property randomly
|
|
set seed based on atom tag
|
|
make atom result independent of what proc owns it
|
|
------------------------------------------------------------------------- */
|
|
|
|
void Set::setrandom(int keyword)
|
|
{
|
|
int i;
|
|
|
|
selection(atom->nlocal);
|
|
RanPark *random = new RanPark(lmp,1);
|
|
double **x = atom->x;
|
|
int seed = ivalue;
|
|
|
|
// set fraction of atom types to newtype
|
|
|
|
if (keyword == TYPE_FRACTION) {
|
|
int nlocal = atom->nlocal;
|
|
|
|
for (i = 0; i < nlocal; i++)
|
|
if (select[i]) {
|
|
random->reset(seed,x[i]);
|
|
if (random->uniform() > fraction) continue;
|
|
atom->type[i] = newtype;
|
|
count++;
|
|
}
|
|
|
|
// set dipole moments to random orientations in 3d or 2d
|
|
// dipole length is determined by dipole type array
|
|
|
|
} else if (keyword == DIPOLE_RANDOM) {
|
|
double **mu = atom->mu;
|
|
int nlocal = atom->nlocal;
|
|
|
|
double msq,scale;
|
|
|
|
if (domain->dimension == 3) {
|
|
for (i = 0; i < nlocal; i++)
|
|
if (select[i]) {
|
|
random->reset(seed,x[i]);
|
|
mu[i][0] = random->uniform() - 0.5;
|
|
mu[i][1] = random->uniform() - 0.5;
|
|
mu[i][2] = random->uniform() - 0.5;
|
|
msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2];
|
|
scale = dvalue/sqrt(msq);
|
|
mu[i][0] *= scale;
|
|
mu[i][1] *= scale;
|
|
mu[i][2] *= scale;
|
|
mu[i][3] = dvalue;
|
|
count++;
|
|
}
|
|
|
|
} else {
|
|
for (i = 0; i < nlocal; i++)
|
|
if (select[i]) {
|
|
random->reset(seed,x[i]);
|
|
mu[i][0] = random->uniform() - 0.5;
|
|
mu[i][1] = random->uniform() - 0.5;
|
|
mu[i][2] = 0.0;
|
|
msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1];
|
|
scale = dvalue/sqrt(msq);
|
|
mu[i][0] *= scale;
|
|
mu[i][1] *= scale;
|
|
mu[i][3] = dvalue;
|
|
count++;
|
|
}
|
|
}
|
|
|
|
// set quaternions to random orientations in 3d or 2d
|
|
// no need to normalize quats since creations algorithms already do
|
|
|
|
} else if (keyword == QUAT_RANDOM) {
|
|
AtomVecEllipsoid *avec_ellipsoid =
|
|
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
|
AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
|
|
int *ellipsoid = atom->ellipsoid;
|
|
int nlocal = atom->nlocal;
|
|
double *quat;
|
|
|
|
if (domain->dimension == 3) {
|
|
double s,t1,t2,theta1,theta2;
|
|
for (i = 0; i < nlocal; i++)
|
|
if (select[i]) {
|
|
if (ellipsoid[i] < 0)
|
|
error->one(FLERR,"Cannot set quaternion for atom "
|
|
"that is not an ellipsoid");
|
|
quat = bonus[ellipsoid[i]].quat;
|
|
random->reset(seed,x[i]);
|
|
s = random->uniform();
|
|
t1 = sqrt(1.0-s);
|
|
t2 = sqrt(s);
|
|
theta1 = 2.0*PI*random->uniform();
|
|
theta2 = 2.0*PI*random->uniform();
|
|
quat[0] = cos(theta2)*t2;
|
|
quat[1] = sin(theta1)*t1;
|
|
quat[2] = cos(theta1)*t1;
|
|
quat[3] = sin(theta2)*t2;
|
|
count++;
|
|
}
|
|
|
|
} else {
|
|
double theta2;
|
|
for (i = 0; i < nlocal; i++)
|
|
if (select[i]) {
|
|
if (ellipsoid[i] < 0)
|
|
error->one(FLERR,"Cannot set quaternion for atom "
|
|
"that is not an ellipsoid");
|
|
quat = bonus[ellipsoid[i]].quat;
|
|
random->reset(seed,x[i]);
|
|
theta2 = PI*random->uniform();
|
|
quat[0] = cos(theta2);
|
|
quat[1] = 0.0;
|
|
quat[2] = 0.0;
|
|
quat[3] = sin(theta2);
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
|
|
delete random;
|
|
}
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
void Set::topology(int keyword)
|
|
{
|
|
int m,atom1,atom2,atom3,atom4;
|
|
|
|
// border swap to acquire ghost atom info
|
|
// enforce PBC before in case atoms are outside box
|
|
// init entire system since comm->exchange is done
|
|
// comm::init needs neighbor::init needs pair::init needs kspace::init, etc
|
|
|
|
if (comm->me == 0 && screen) fprintf(screen," system init for set ...\n");
|
|
lmp->init();
|
|
|
|
if (domain->triclinic) domain->x2lamda(atom->nlocal);
|
|
domain->pbc();
|
|
domain->reset_box();
|
|
comm->setup();
|
|
comm->exchange();
|
|
comm->borders();
|
|
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
|
|
|
|
// select both owned and ghost atoms
|
|
|
|
selection(atom->nlocal + atom->nghost);
|
|
|
|
// for BOND, each of 2 atoms must be in group
|
|
|
|
if (keyword == BOND) {
|
|
int nlocal = atom->nlocal;
|
|
for (int i = 0; i < nlocal; i++)
|
|
for (m = 0; m < atom->num_bond[i]; m++) {
|
|
atom1 = atom->map(atom->bond_atom[i][m]);
|
|
if (atom1 == -1) error->one(FLERR,"Bond atom missing in set command");
|
|
if (select[i] && select[atom1]) {
|
|
atom->bond_type[i][m] = ivalue;
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// for ANGLE, each of 3 atoms must be in group
|
|
|
|
if (keyword == ANGLE) {
|
|
int nlocal = atom->nlocal;
|
|
for (int i = 0; i < nlocal; i++)
|
|
for (m = 0; m < atom->num_angle[i]; m++) {
|
|
atom1 = atom->map(atom->angle_atom1[i][m]);
|
|
atom2 = atom->map(atom->angle_atom2[i][m]);
|
|
atom3 = atom->map(atom->angle_atom3[i][m]);
|
|
if (atom1 == -1 || atom2 == -1 || atom3 == -1)
|
|
error->one(FLERR,"Angle atom missing in set command");
|
|
if (select[atom1] && select[atom2] && select[atom3]) {
|
|
atom->angle_type[i][m] = ivalue;
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// for DIHEDRAL, each of 4 atoms must be in group
|
|
|
|
if (keyword == DIHEDRAL) {
|
|
int nlocal = atom->nlocal;
|
|
for (int i = 0; i < nlocal; i++)
|
|
for (m = 0; m < atom->num_dihedral[i]; m++) {
|
|
atom1 = atom->map(atom->dihedral_atom1[i][m]);
|
|
atom2 = atom->map(atom->dihedral_atom2[i][m]);
|
|
atom3 = atom->map(atom->dihedral_atom3[i][m]);
|
|
atom4 = atom->map(atom->dihedral_atom4[i][m]);
|
|
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
|
|
error->one(FLERR,"Dihedral atom missing in set command");
|
|
if (select[atom1] && select[atom2] && select[atom3] && select[atom4]) {
|
|
atom->dihedral_type[i][m] = ivalue;
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// for IMPROPER, each of 4 atoms must be in group
|
|
|
|
if (keyword == IMPROPER) {
|
|
int nlocal = atom->nlocal;
|
|
for (int i = 0; i < nlocal; i++)
|
|
for (m = 0; m < atom->num_improper[i]; m++) {
|
|
atom1 = atom->map(atom->improper_atom1[i][m]);
|
|
atom2 = atom->map(atom->improper_atom2[i][m]);
|
|
atom3 = atom->map(atom->improper_atom3[i][m]);
|
|
atom4 = atom->map(atom->improper_atom4[i][m]);
|
|
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
|
|
error->one(FLERR,"Improper atom missing in set command");
|
|
if (select[atom1] && select[atom2] && select[atom3] && select[atom4]) {
|
|
atom->improper_type[i][m] = ivalue;
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
}
|