ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14506 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2016-01-27 20:33:55 +00:00
parent dca90d44b7
commit fa7543b714
7 changed files with 150 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
src/GRANULAR/fix_pour.cpp
View File

@ -156,7 +156,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
memory->create(coords,natom_max,4,"pour:coords");
memory->create(imageflags,natom_max,"pour:imageflags");
// find current max atom and molecule IDs if necessary
// find max atom and molecule IDs just once
if (idnext) find_maxid();
@ -385,7 +385,7 @@ void FixPour::pre_exchange()
atom->nghost = 0;
atom->avec->clear_bonus();
// find current max atom and molecule IDs if necessary
// find current max atom and molecule IDs on every insertion step
if (!idnext) find_maxid();
@ -737,12 +737,19 @@ void FixPour::find_maxid()
return 1 if yes, 0 if no
for ATOM mode, use delta with maximum size for inserted atoms
for MOLECULE mode, use delta with max radius of inserted molecules
if ignore line/tri set, ignore line or tri particles
account for PBC in overlap decision via outside() and minimum_image()
------------------------------------------------------------------------- */
int FixPour::overlap(int i)
{
double delta;
if (ignoreflag) {
if (ignoreline && atom->line[i] >= 0) return 0;
if (ignoretri && atom->tri[i] >= 0) return 0;
}
if (mode == ATOM) delta = atom->radius[i] + radius_max;
else delta = atom->radius[i] + molradius_max;
@ -860,6 +867,7 @@ void FixPour::options(int narg, char **arg)
shakeflag = 0;
idshake = NULL;
idnext = 0;
ignoreflag = ignoreline = ignoretri = 0;
dstyle = ONE;
radius_max = radius_one = 0.5;
radius_poly = frac_poly = NULL;
@ -925,6 +933,13 @@ void FixPour::options(int narg, char **arg)
else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1;
else error->all(FLERR,"Illegal fix pour command");
iarg += 2;
} else if (strcmp(arg[iarg],"ignore") == 0) {
if (atom->line_flag) ignoreline = 1;
if (atom->tri_flag) ignoretri = 1;
if (ignoreline || ignoretri) ignoreflag = 1;
iarg += 1;
} else if (strcmp(arg[iarg],"diam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
if (strcmp(arg[iarg+1],"one") == 0) {
@ -938,6 +953,7 @@ void FixPour::options(int narg, char **arg)
dstyle = RANGE;
radius_lo = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_hi = 0.5 * force->numeric(FLERR,arg[iarg+3]);
if (radius_lo > radius_hi) error->all(FLERR,"Illegal fix pour command");
radius_max = radius_hi;
iarg += 4;
} else if (strcmp(arg[iarg+1],"poly") == 0) {
@ -968,6 +984,7 @@ void FixPour::options(int narg, char **arg)
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
density_lo = force->numeric(FLERR,arg[iarg+1]);
density_hi = force->numeric(FLERR,arg[iarg+2]);
if (density_lo > density_hi) error->all(FLERR,"Illegal fix pour command");
iarg += 3;
} else if (strcmp(arg[iarg],"vol") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
@ -985,6 +1002,8 @@ void FixPour::options(int narg, char **arg)
vxhi = force->numeric(FLERR,arg[iarg+2]);
vylo = force->numeric(FLERR,arg[iarg+3]);
vyhi = force->numeric(FLERR,arg[iarg+4]);
if (vxlo > vxhi || vylo > vyhi)
error->all(FLERR,"Illegal fix pour command");
vz = force->numeric(FLERR,arg[iarg+5]);
iarg += 6;
} else {
@ -993,6 +1012,7 @@ void FixPour::options(int narg, char **arg)
vxhi = force->numeric(FLERR,arg[iarg+2]);
vy = force->numeric(FLERR,arg[iarg+3]);
vz = 0.0;
if (vxlo > vxhi) error->all(FLERR,"Illegal fix pour command");
iarg += 4;
}
} else error->all(FLERR,"Illegal fix pour command");

112
src/atom_vec_tri.cpp
View File

@ -41,7 +41,7 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
size_forward = 7;
size_reverse = 6;
size_border = 26;
size_velocity = 6;
size_velocity = 9;
size_data_atom = 8;
size_data_vel = 7;
size_data_bonus = 10;
@ -49,7 +49,9 @@ AtomVecTri::AtomVecTri(LAMMPS *lmp) : AtomVec(lmp)
atom->tri_flag = 1;
atom->molecule_flag = atom->rmass_flag = 1;
atom->radius_flag = atom->angmom_flag = atom->torque_flag = 1;
atom->radius_flag = atom->omega_flag = atom->angmom_flag = 1;
atom->torque_flag = 1;
atom->sphere_flag = 1;
nlocal_bonus = nghost_bonus = nmax_bonus = 0;
bonus = NULL;
@ -100,6 +102,7 @@ void AtomVecTri::grow(int n)
molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
radius = memory->grow(atom->radius,nmax,"atom:radius");
omega = memory->grow(atom->omega,nmax,3,"atom:omega");
angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
tri = memory->grow(atom->tri,nmax,"atom:tri");
@ -119,7 +122,8 @@ void AtomVecTri::grow_reset()
mask = atom->mask; image = atom->image;
x = atom->x; v = atom->v; f = atom->f;
molecule = atom->molecule; rmass = atom->rmass;
radius = atom->radius; angmom = atom->angmom; torque = atom->torque;
radius = atom->radius; omega = atom->omega;
angmom = atom->angmom; torque = atom->torque;
tri = atom->tri;
}
@ -158,6 +162,9 @@ void AtomVecTri::copy(int i, int j, int delflag)
molecule[j] = molecule[i];
rmass[j] = rmass[i];
radius[j] = radius[i];
omega[j][0] = omega[i][0];
omega[j][1] = omega[i][1];
omega[j][2] = omega[i][2];
angmom[j][0] = angmom[i][0];
angmom[j][1] = angmom[i][1];
angmom[j][2] = angmom[i][2];
@ -343,6 +350,9 @@ int AtomVecTri::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
buf[m++] = angmom[j][0];
buf[m++] = angmom[j][1];
buf[m++] = angmom[j][2];
@ -373,6 +383,9 @@ int AtomVecTri::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
buf[m++] = angmom[j][0];
buf[m++] = angmom[j][1];
buf[m++] = angmom[j][2];
@ -402,6 +415,9 @@ int AtomVecTri::pack_comm_vel(int n, int *list, double *buf,
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
buf[m++] = angmom[j][0];
buf[m++] = angmom[j][1];
buf[m++] = angmom[j][2];
@ -479,6 +495,9 @@ void AtomVecTri::unpack_comm_vel(int n, int first, double *buf)
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
omega[i][0] = buf[m++];
omega[i][1] = buf[m++];
omega[i][2] = buf[m++];
angmom[i][0] = buf[m++];
angmom[i][1] = buf[m++];
angmom[i][2] = buf[m++];
@ -728,6 +747,9 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf,
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
buf[m++] = angmom[j][0];
buf[m++] = angmom[j][1];
buf[m++] = angmom[j][2];
@ -782,6 +804,9 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf,
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
buf[m++] = angmom[j][0];
buf[m++] = angmom[j][1];
buf[m++] = angmom[j][2];
@ -835,6 +860,9 @@ int AtomVecTri::pack_border_vel(int n, int *list, double *buf,
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
buf[m++] = omega[j][0];
buf[m++] = omega[j][1];
buf[m++] = omega[j][2];
buf[m++] = angmom[j][0];
buf[m++] = angmom[j][1];
buf[m++] = angmom[j][2];
@ -1002,6 +1030,9 @@ void AtomVecTri::unpack_border_vel(int n, int first, double *buf)
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
omega[i][0] = buf[m++];
omega[i][1] = buf[m++];
omega[i][2] = buf[m++];
angmom[i][0] = buf[m++];
angmom[i][1] = buf[m++];
angmom[i][2] = buf[m++];
@ -1082,6 +1113,9 @@ int AtomVecTri::pack_exchange(int i, double *buf)
buf[m++] = ubuf(molecule[i]).d;
buf[m++] = rmass[i];
buf[m++] = radius[i];
buf[m++] = omega[i][0];
buf[m++] = omega[i][1];
buf[m++] = omega[i][2];
buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1];
buf[m++] = angmom[i][2];
@ -1143,6 +1177,9 @@ int AtomVecTri::unpack_exchange(double *buf)
molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
rmass[nlocal] = buf[m++];
radius[nlocal] = buf[m++];
omega[nlocal][0] = buf[m++];
omega[nlocal][1] = buf[m++];
omega[nlocal][2] = buf[m++];
angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++];
angmom[nlocal][2] = buf[m++];
@ -1197,8 +1234,8 @@ int AtomVecTri::size_restart()
int n = 0;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (tri[i] >= 0) n += 34;
else n += 18;
if (tri[i] >= 0) n += 37;
else n += 21;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
@ -1231,6 +1268,9 @@ int AtomVecTri::pack_restart(int i, double *buf)
buf[m++] = ubuf(molecule[i]).d;
buf[m++] = rmass[i];
buf[m++] = radius[i];
buf[m++] = omega[i][0];
buf[m++] = omega[i][1];
buf[m++] = omega[i][2];
buf[m++] = angmom[i][0];
buf[m++] = angmom[i][1];
buf[m++] = angmom[i][2];
@ -1298,6 +1338,9 @@ int AtomVecTri::unpack_restart(double *buf)
molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
rmass[nlocal] = buf[m++];
radius[nlocal] = buf[m++];
omega[nlocal][0] = buf[m++];
omega[nlocal][1] = buf[m++];
omega[nlocal][2] = buf[m++];
angmom[nlocal][0] = buf[m++];
angmom[nlocal][1] = buf[m++];
angmom[nlocal][2] = buf[m++];
@ -1366,6 +1409,9 @@ void AtomVecTri::create_atom(int itype, double *coord)
molecule[nlocal] = 0;
radius[nlocal] = 0.5;
rmass[nlocal] = 4.0*MY_PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal];
omega[nlocal][0] = 0.0;
omega[nlocal][1] = 0.0;
omega[nlocal][2] = 0.0;
angmom[nlocal][0] = 0.0;
angmom[nlocal][1] = 0.0;
angmom[nlocal][2] = 0.0;
@ -1415,6 +1461,9 @@ void AtomVecTri::data_atom(double *coord, imageint imagetmp, char **values)
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
v[nlocal][2] = 0.0;
omega[nlocal][0] = 0.0;
omega[nlocal][1] = 0.0;
omega[nlocal][2] = 0.0;
angmom[nlocal][0] = 0.0;
angmom[nlocal][1] = 0.0;
angmom[nlocal][2] = 0.0;
@ -1591,21 +1640,27 @@ void AtomVecTri::data_vel(int m, char **values)
v[m][0] = atof(values[0]);
v[m][1] = atof(values[1]);
v[m][2] = atof(values[2]);
angmom[m][0] = atof(values[3]);
angmom[m][1] = atof(values[4]);
angmom[m][2] = atof(values[5]);
omega[m][0] = atof(values[3]);
omega[m][1] = atof(values[4]);
omega[m][2] = atof(values[5]);
angmom[m][0] = atof(values[6]);
angmom[m][1] = atof(values[7]);
angmom[m][2] = atof(values[8]);
}
/* ----------------------------------------------------------------------
unpack hybrid quantities from one tri in Velocities section of data file
unpack hybrid quantities from one line in Velocities section of data file
------------------------------------------------------------------------- */
int AtomVecTri::data_vel_hybrid(int m, char **values)
{
angmom[m][0] = atof(values[0]);
angmom[m][1] = atof(values[1]);
angmom[m][2] = atof(values[2]);
return 3;
omega[m][0] = atof(values[0]);
omega[m][1] = atof(values[1]);
omega[m][2] = atof(values[2]);
angmom[m][0] = atof(values[3]);
angmom[m][1] = atof(values[4]);
angmom[m][2] = atof(values[5]);
return 6;
}
/* ----------------------------------------------------------------------
@ -1703,9 +1758,12 @@ void AtomVecTri::pack_vel(double **buf)
buf[i][1] = v[i][0];
buf[i][2] = v[i][1];
buf[i][3] = v[i][2];
buf[i][4] = angmom[i][0];
buf[i][5] = angmom[i][1];
buf[i][6] = angmom[i][2];
buf[i][4] = omega[i][0];
buf[i][5] = omega[i][1];
buf[i][6] = omega[i][2];
buf[i][7] = angmom[i][0];
buf[i][8] = angmom[i][1];
buf[i][9] = angmom[i][2];
}
}
@ -1715,10 +1773,13 @@ void AtomVecTri::pack_vel(double **buf)
int AtomVecTri::pack_vel_hybrid(int i, double *buf)
{
buf[0] = angmom[i][0];
buf[1] = angmom[i][1];
buf[2] = angmom[i][2];
return 3;
buf[0] = omega[i][0];
buf[1] = omega[i][1];
buf[2] = omega[i][2];
buf[3] = angmom[i][0];
buf[4] = angmom[i][1];
buf[5] = angmom[i][2];
return 6;
}
/* ----------------------------------------------------------------------
@ -1729,9 +1790,10 @@ void AtomVecTri::write_vel(FILE *fp, int n, double **buf)
{
for (int i = 0; i < n; i++)
fprintf(fp,TAGINT_FORMAT
" %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n",
" %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e "
"%-1.16e %-1.16e %-1.16e\n",
(tagint) ubuf(buf[i][0]).i,buf[i][1],buf[i][2],buf[i][3],
buf[i][4],buf[i][5],buf[i][6]);
buf[i][4],buf[i][5],buf[i][6],buf[i][7],buf[i][8],buf[i][9]);
}
/* ----------------------------------------------------------------------
@ -1740,8 +1802,9 @@ void AtomVecTri::write_vel(FILE *fp, int n, double **buf)
int AtomVecTri::write_vel_hybrid(FILE *fp, double *buf)
{
fprintf(fp," %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2]);
return 3;
fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
buf[0],buf[1],buf[2],buf[3],buf[4],buf[5]);
return 6;
}
/* ----------------------------------------------------------------------
@ -1763,6 +1826,7 @@ bigint AtomVecTri::memory_usage()
if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
if (atom->memcheck("radius")) bytes += memory->usage(radius,nmax);
if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
if (atom->memcheck("torque")) bytes +=
memory->usage(torque,nmax*comm->nthreads,3);

2
src/atom_vec_tri.h
View File

@ -92,7 +92,7 @@ class AtomVecTri : public AtomVec {
double **x,**v,**f;
tagint *molecule;
double *rmass,*radius;
double **angmom,**torque;
double **omega,**angmom,**torque;
int *tri;
int nlocal_bonus,nghost_bonus,nmax_bonus;

14
src/displace_atoms.cpp
View File

@ -232,8 +232,16 @@ void DisplaceAtoms::command(int narg, char **arg)
runit[2] = axis[2]/len;
double angle = MY_PI*theta/180.0;
double sine = sin(angle);
double cosine = cos(angle);
double sine = sin(angle);
double qcosine = cos(0.5*angle);
double qsine = sin(0.5*angle);
qrotate[0] = qcosine;
qrotate[1] = runit[0]*qsine;
qrotate[2] = runit[1]*qsine;
qrotate[3] = runit[2]*qsine;
double ddotr;
// flags for additional orientation info stored by some atom styles
@ -304,10 +312,6 @@ void DisplaceAtoms::command(int narg, char **arg)
else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat;
if (quat) {
qrotate[0] = cosine;
qrotate[1] = runit[0]*sine;
qrotate[2] = runit[1]*sine;
qrotate[3] = runit[2]*sine;
MathExtra::quatquat(qrotate,quat,qnew);
quat[0] = qnew[0];
quat[1] = qnew[1];

20
src/dump_image.cpp
View File

@ -178,11 +178,13 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
iarg += 3;
} else if (strcmp(arg[iarg],"tri") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command");
if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
triflag = YES;
if (strcmp(arg[iarg+1],"type") == 0) tcolor = TYPE;
else error->all(FLERR,"Illegal dump image command");
iarg += 2;
tstyle = force->inumeric(FLERR,arg[iarg+2]);
tdiamvalue = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"body") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command");
@ -806,12 +808,19 @@ void DumpImage::create_image()
}
// render atoms that are triangles
// tstyle = 1 for tri only, 2 for edges only, 3 for both
if (triflag) {
int tridraw = 1;
if (tstyle == 2) tridraw = 0;
int edgedraw = 1;
if (tstyle == 1) edgedraw = 0;
double **x = atom->x;
int *tri = atom->tri;
int *type = atom->type;
for (i = 0; i < nchoose; i++) {
j = clist[i];
if (tri[j] < 0) continue;
@ -828,7 +837,12 @@ void DumpImage::create_image()
MathExtra::add3(pt2,x[i],pt2);
MathExtra::add3(pt3,x[i],pt3);
image->draw_triangle(pt1,pt2,pt3,color);
if (tridraw) image->draw_triangle(pt1,pt2,pt3,color);
if (edgedraw) {
image->draw_cylinder(pt1,pt2,color,tdiamvalue,3);
image->draw_cylinder(pt2,pt3,color,tdiamvalue,3);
image->draw_cylinder(pt3,pt1,color,tdiamvalue,3);
}
}
}

3
src/dump_image.h
View File

@ -45,7 +45,8 @@ class DumpImage : public DumpCustom {
int lcolor,ldiam; // what determines color/diam of lines
double ldiamvalue; // line diameter value
int triflag; // 0/1 for draw atoms as triangles
int tcolor; // what determines color of tris
int tcolor,tstyle; // what determines color/style of tris
double tdiamvalue; // tri edge diameter value
int bodyflag; // 0/1 for draw atoms as bodies
int bodycolor; // what determines color of bodies
double bodyflag1,bodyflag2; // user-specified params for drawing bodies

25
src/fix_move.cpp
View File

@ -193,11 +193,6 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
"Fix move cannot define z or vz variable for 2d problem");
}
if (atom->angmom_flag && comm->me == 0)
error->warning(FLERR,"Fix move does not update angular momentum");
if (atom->ellipsoid_flag && comm->me == 0)
error->warning(FLERR,"Fix move does not update quaternions");
// setup scaling and apply scaling factors to velocity & amplitude
if ((mstyle == LINEAR || mstyle == WIGGLE || mstyle == ROTATE) &&
@ -609,8 +604,15 @@ void FixMove::initial_integrate(int vflag)
} else if (mstyle == ROTATE) {
double arg = omega_rotate * delta;
double sine = sin(arg);
double cosine = cos(arg);
double sine = sin(arg);
double qcosine = cos(0.5*arg);
double qsine = sin(0.5*arg);
qrotate[0] = qcosine;
qrotate[1] = runit[0]*qsine;
qrotate[2] = runit[1]*qsine;
qrotate[3] = runit[2]*qsine;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
@ -646,12 +648,13 @@ void FixMove::initial_integrate(int vflag)
if (extra_flag) {
// omega for spheres and lines
// omega for spheres, lines, tris
if (omega_flag) {
flag = 0;
if (radius_flag && radius[i] > 0.0) flag = 1;
if (line_flag && line[i] >= 0.0) flag = 1;
if (tri_flag && tri[i] >= 0.0) flag = 1;
if (flag) {
omega[i][0] = omega_rotate*runit[0];
omega[i][1] = omega_rotate*runit[1];
@ -706,13 +709,7 @@ void FixMove::initial_integrate(int vflag)
quat = avec_tri->bonus[tri[i]].quat;
else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat;
if (quat) {
qrotate[0] = cosine;
qrotate[1] = runit[0]*sine;
qrotate[2] = runit[1]*sine;
qrotate[3] = runit[2]*sine;
MathExtra::quatquat(qrotate,qoriginal[i],quat);
}
if (quat) MathExtra::quatquat(qrotate,qoriginal[i],quat);
}
}