From 09fb68df71adf6b39feaa3dfb8c0be9bb53d674e Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 14 Jan 2016 22:33:46 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14440 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/pair_line_lj.cpp | 93 +++++++----- src/ASPHERE/pair_line_lj.h | 8 +- src/ASPHERE/pair_tri_lj.cpp | 2 +- src/displace_atoms.cpp | 70 ++++++++- src/dump_image.cpp | 200 +++++++++++++++++++----- src/dump_image.h | 20 ++- src/fix_move.cpp | 286 ++++++++++++++++++++++++++++++++--- src/fix_move.h | 12 +- src/math_extra.cpp | 1 + src/math_extra.h | 13 +- 10 files changed, 588 insertions(+), 117 deletions(-) diff --git a/src/ASPHERE/pair_line_lj.cpp b/src/ASPHERE/pair_line_lj.cpp index 2124b455f6..737aefa18d 100644 --- a/src/ASPHERE/pair_line_lj.cpp +++ b/src/ASPHERE/pair_line_lj.cpp @@ -52,7 +52,10 @@ PairLineLJ::~PairLineLJ() memory->destroy(setflag); memory->destroy(cutsq); + memory->destroy(subsize); memory->destroy(cut); + memory->destroy(cutsub); + memory->destroy(cutsubsq); memory->destroy(epsilon); memory->destroy(sigma); memory->destroy(lj1); @@ -66,7 +69,7 @@ PairLineLJ::~PairLineLJ() void PairLineLJ::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype,tmp; int ni,nj,npi,npj,ifirst,jfirst; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,term1,term2,sig,sig3,forcelj; @@ -130,10 +133,10 @@ void PairLineLJ::compute(int eflag, int vflag) evdwl = 0.0; if (line[i] >= 0 && line[j] >= 0) { - if (dnum[i] == 0) discretize(i,sigma[itype][itype]); + if (dnum[i] == 0) discretize(i,subsize[itype]); npi = dnum[i]; ifirst = dfirst[i]; - if (dnum[j] == 0) discretize(j,sigma[jtype][jtype]); + if (dnum[j] == 0) discretize(j,subsize[jtype]); npj = dnum[j]; jfirst = dfirst[j]; @@ -154,7 +157,11 @@ void PairLineLJ::compute(int eflag, int vflag) dely = xi[1] - xj[1]; rsq = delx*delx + dely*dely; - sig = 0.5 * (discrete[ifirst+ni].sigma+discrete[jfirst+nj].sigma); + // skip this pair of sub-particles if outside sub cutoff + + if (rsq >= cutsubsq[itype][jtype]) continue; + + sig = sigma[itype][jtype]; sig3 = sig*sig*sig; term2 = 24.0*epsilon[itype][jtype] * sig3*sig3; term1 = 2.0 * term2 * sig3*sig3; @@ -167,16 +174,15 @@ void PairLineLJ::compute(int eflag, int vflag) fi[0] = delx*fpair; fi[1] = dely*fpair; + f[i][0] += fi[0]; f[i][1] += fi[1]; torque[i][2] += dxi*fi[1] - dyi*fi[0]; if (newton_pair || j < nlocal) { - fj[0] = -delx*fpair; - fj[1] = -dely*fpair; - f[j][0] += fj[0]; - f[j][1] += fj[1]; - torque[j][2] += dxj*fj[1] - dyj*fj[0]; + f[j][0] -= fi[0]; + f[j][1] -= fi[1]; + torque[j][2] -= dxj*fi[1] - dyj*fi[0]; } } } @@ -185,7 +191,7 @@ void PairLineLJ::compute(int eflag, int vflag) // convert line into Np particles based on sigma and line length } else if (line[i] >= 0) { - if (dnum[i] == 0) discretize(i,sigma[itype][itype]); + if (dnum[i] == 0) discretize(i,subsize[itype]); npi = dnum[i]; ifirst = dfirst[i]; @@ -202,7 +208,11 @@ void PairLineLJ::compute(int eflag, int vflag) dely = xi[1] - xj[1]; rsq = delx*delx + dely*dely; - sig = 0.5 * (discrete[ifirst+ni].sigma+sigma[jtype][jtype]); + // skip this pair of sub-particles if outside sub cutoff + + if (rsq >= cutsubsq[itype][jtype]) continue; + + sig = sigma[itype][jtype]; sig3 = sig*sig*sig; term2 = 24.0*epsilon[itype][jtype] * sig3*sig3; term1 = 2.0 * term2 * sig3*sig3; @@ -220,10 +230,8 @@ void PairLineLJ::compute(int eflag, int vflag) torque[i][2] += dxi*fi[1] - dyi*fi[0]; if (newton_pair || j < nlocal) { - fj[0] = -delx*fpair; - fj[1] = -dely*fpair; - f[j][0] += fj[0]; - f[j][1] += fj[1]; + f[j][0] -= fi[0]; + f[j][1] -= fi[1]; } } @@ -231,7 +239,7 @@ void PairLineLJ::compute(int eflag, int vflag) // convert line into Np particles based on sigma and line length } else if (line[j] >= 0) { - if (dnum[j] == 0) discretize(j,sigma[jtype][jtype]); + if (dnum[j] == 0) discretize(j,subsize[jtype]); npj = dnum[j]; jfirst = dfirst[j]; @@ -248,7 +256,11 @@ void PairLineLJ::compute(int eflag, int vflag) dely = xi[1] - xj[1]; rsq = delx*delx + dely*dely; - sig = 0.5 * (sigma[itype][itype]+discrete[jfirst+nj].sigma); + // skip this pair of sub-particles if outside sub cutoff + + if (rsq >= cutsubsq[itype][jtype]) continue; + + sig = sigma[itype][jtype]; sig3 = sig*sig*sig; term2 = 24.0*epsilon[itype][jtype] * sig3*sig3; term1 = 2.0 * term2 * sig3*sig3; @@ -265,11 +277,9 @@ void PairLineLJ::compute(int eflag, int vflag) f[i][1] += fi[1]; if (newton_pair || j < nlocal) { - f[j][0] += fj[0]; - f[j][1] += fj[1]; - fj[0] = -delx*fpair; - fj[1] = -dely*fpair; - torque[j][2] += dxj*fj[1] - dyj*fj[0]; + f[j][0] -= fi[0]; + f[j][1] -= fi[1]; + torque[j][2] -= dxj*fi[1] - dyj*fi[0]; } } @@ -318,7 +328,10 @@ void PairLineLJ::allocate() memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(subsize,n+1,"pair:subsize"); memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(cutsub,n+1,n+1,"pair:cutsub"); + memory->create(cutsubsq,n+1,n+1,"pair:cutsubsq"); memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(sigma,n+1,n+1,"pair:sigma"); memory->create(lj1,n+1,n+1,"pair:lj1"); @@ -353,7 +366,7 @@ void PairLineLJ::settings(int narg, char **arg) void PairLineLJ::coeff(int narg, char **arg) { - if (narg < 4 || narg > 5) + if (narg < 7 || narg > 8) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -361,17 +374,23 @@ void PairLineLJ::coeff(int narg, char **arg) force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); - double epsilon_one = force->numeric(FLERR,arg[2]); - double sigma_one = force->numeric(FLERR,arg[3]); + double size_itype = force->numeric(FLERR,arg[2]); + double size_jtype = force->numeric(FLERR,arg[3]); + double epsilon_one = force->numeric(FLERR,arg[4]); + double sigma_one = force->numeric(FLERR,arg[5]); + double cutsub_one = force->numeric(FLERR,arg[6]); double cut_one = cut_global; - if (narg == 5) cut_one = force->numeric(FLERR,arg[4]); + if (narg == 8) cut_one = force->numeric(FLERR,arg[7]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { + subsize[i] = size_itype; + subsize[j] = size_jtype; epsilon[i][j] = epsilon_one; sigma[i][j] = sigma_one; + cutsub[i][j] = cutsub_one; cut[i][j] = cut_one; setflag[i][j] = 1; count++; @@ -399,12 +418,9 @@ void PairLineLJ::init_style() double PairLineLJ::init_one(int i, int j) { - if (setflag[i][j] == 0) { - epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], - sigma[i][i],sigma[j][j]); - sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); - } + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + cutsubsq[i][j] = cutsub[i][j] * cutsub[i][j]; lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); @@ -413,6 +429,8 @@ double PairLineLJ::init_one(int i, int j) epsilon[j][i] = epsilon[i][j]; sigma[j][i] = sigma[i][j]; + cutsubsq[j][i] = cutsubsq[i][j]; + lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; @@ -422,16 +440,17 @@ double PairLineLJ::init_one(int i, int j) } /* ---------------------------------------------------------------------- - discretize line segment I into N sub-segments no more than sigma in length - store new discrete particles in Discrete list + discretize line segment I into N sub-particles with <= size separation + store displacement dx,dy of discrete particles in Discrete list ------------------------------------------------------------------------- */ -void PairLineLJ::discretize(int i, double sigma) +void PairLineLJ::discretize(int i, double size) { AtomVecLine::Bonus *bonus = avec->bonus; double length = bonus[atom->line[i]].length; double theta = bonus[atom->line[i]].theta; - int n = static_cast (length/sigma) + 1; + int n = static_cast (length/size) + 1; + dnum[i] = n; dfirst[i] = ndiscrete; @@ -441,14 +460,12 @@ void PairLineLJ::discretize(int i, double sigma) memory->srealloc(discrete,dmax*sizeof(Discrete),"pair:discrete"); } - sigma = length/n; double delta; for (int m = 0; m < n; m++) { delta = -0.5 + (2*m+1)/(2.0*n); discrete[ndiscrete].dx = delta*length*cos(theta); discrete[ndiscrete].dy = delta*length*sin(theta); - discrete[ndiscrete].sigma = sigma; ndiscrete++; } } diff --git a/src/ASPHERE/pair_line_lj.h b/src/ASPHERE/pair_line_lj.h index 7d99e836f0..0eb5bb0e09 100644 --- a/src/ASPHERE/pair_line_lj.h +++ b/src/ASPHERE/pair_line_lj.h @@ -36,14 +36,16 @@ class PairLineLJ : public Pair { protected: double cut_global; + double *subsize; + double **epsilon,**sigma,**cutsub,**cutsubsq; double **cut; - double **epsilon,**sigma; - double **lj1,**lj2,**lj3,**lj4; + double **lj1,**lj2,**lj3,**lj4; // for sphere/sphere interactions class AtomVecLine *avec; + double *size; // per-type size of sub-particles to tile line segment + struct Discrete { double dx,dy; - double sigma; }; Discrete *discrete; // list of all discretes for all lines int ndiscrete; // number of discretes in list diff --git a/src/ASPHERE/pair_tri_lj.cpp b/src/ASPHERE/pair_tri_lj.cpp index 4823664de8..ef5123dc11 100644 --- a/src/ASPHERE/pair_tri_lj.cpp +++ b/src/ASPHERE/pair_tri_lj.cpp @@ -513,7 +513,7 @@ double PairTriLJ::init_one(int i, int j) ------------------------------------------------------------------------- */ void PairTriLJ::discretize(int i, double sigma, - double *c1, double *c2, double *c3) + double *c1, double *c2, double *c3) { double centroid[3],dc1[3],dc2[3],dc3[3]; diff --git a/src/displace_atoms.cpp b/src/displace_atoms.cpp index 3c86702a9d..8d2e501641 100644 --- a/src/displace_atoms.cpp +++ b/src/displace_atoms.cpp @@ -24,11 +24,16 @@ #include "group.h" #include "math_const.h" #include "random_park.h" -#include "error.h" #include "force.h" #include "input.h" #include "variable.h" +#include "atom_vec_ellipsoid.h" +#include "atom_vec_line.h" +#include "atom_vec_tri.h" +#include "atom_vec_body.h" +#include "math_extra.h" #include "memory.h" +#include "error.h" using namespace LAMMPS_NS; using namespace MathConst; @@ -203,8 +208,10 @@ void DisplaceAtoms::command(int narg, char **arg) // X = P + C + A cos(theta) + B sin(theta) if (style == ROTATE) { - double axis[3],point[3]; + double theta_new; + double axis[3],point[3],qrotate[4],qnew[4]; double a[3],b[3],c[3],d[3],disp[3],runit[3]; + double *quat; int dim = domain->dimension; point[0] = xscale*force->numeric(FLERR,arg[2]); @@ -224,11 +231,36 @@ void DisplaceAtoms::command(int narg, char **arg) runit[1] = axis[1]/len; runit[2] = axis[2]/len; - double sine = sin(MY_PI*theta/180.0); - double cosine = cos(MY_PI*theta/180.0); + double angle = MY_PI*theta/180.0; + double sine = sin(angle); + double cosine = cos(angle); double ddotr; + // flags for additional orientation info stored by some atom styles + + int ellipsoid_flag = atom->ellipsoid_flag; + int line_flag = atom->line_flag; + int tri_flag = atom->tri_flag; + int body_flag = atom->body_flag; + + int theta_flag = 0; + int quat_flag = 0; + if (line_flag) theta_flag = 1; + if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1; + + // AtomVec pointers to retrieve per-atom storage of extra quantities + + AtomVecEllipsoid *avec_ellipsoid = + (AtomVecEllipsoid *) atom->style_match("ellipsoid"); + AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line"); + AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri"); + AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body"); + double **x = atom->x; + int *ellipsoid = atom->ellipsoid; + int *line = atom->line; + int *tri = atom->tri; + int *body = atom->body; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -253,6 +285,36 @@ void DisplaceAtoms::command(int narg, char **arg) x[i][0] = point[0] + c[0] + disp[0]; x[i][1] = point[1] + c[1] + disp[1]; if (dim == 3) x[i][2] = point[2] + c[2] + disp[2]; + + // theta for lines + + if (theta_flag && line[i] >= 0.0) { + theta_new = fmod(avec_line->bonus[line[i]].theta+angle,MY_2PI); + avec_line->bonus[atom->line[i]].theta = theta_new; + } + + // quats for ellipsoids, tris, and bodies + + if (quat_flag) { + quat = NULL; + if (ellipsoid_flag && ellipsoid[i] >= 0) + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (tri_flag && tri[i] >= 0) + 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,quat,qnew); + quat[0] = qnew[0]; + quat[1] = qnew[1]; + quat[2] = qnew[2]; + quat[3] = qnew[3]; + } + } } } } diff --git a/src/dump_image.cpp b/src/dump_image.cpp index fa0fa12be4..e587e661cb 100644 --- a/src/dump_image.cpp +++ b/src/dump_image.cpp @@ -18,7 +18,8 @@ #include "dump_image.h" #include "image.h" #include "atom.h" -#include "atom_vec.h" +#include "atom_vec_line.h" +#include "atom_vec_tri.h" #include "atom_vec_body.h" #include "body.h" #include "molecule.h" @@ -29,6 +30,7 @@ #include "input.h" #include "variable.h" #include "math_const.h" +#include "math_extra.h" #include "error.h" #include "memory.h" @@ -106,7 +108,7 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : // set defaults for optional args atomflag = YES; - bodyflag = NO; + lineflag = triflag = bodyflag = NO; if (atom->nbondtypes == 0) bondflag = NO; else { bondflag = YES; @@ -132,28 +134,19 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : int iarg = ioptional; while (iarg < narg) { - if (strcmp(arg[iarg],"adiam") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); - adiam = NUMERIC; - adiamvalue = force->numeric(FLERR,arg[iarg+1]); - if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command"); - iarg += 2; - - } else if (strcmp(arg[iarg],"atom") == 0) { + if (strcmp(arg[iarg],"atom") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); if (strcmp(arg[iarg+1],"yes") == 0) atomflag = YES; else if (strcmp(arg[iarg+1],"no") == 0) atomflag = NO; else error->all(FLERR,"Illegal dump image command"); iarg += 2; - } else if (strcmp(arg[iarg],"body") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command"); - if (strcmp(arg[iarg+1],"yes") == 0) bodyflag = YES; - else if (strcmp(arg[iarg+1],"no") == 0) bodyflag = NO; - else error->all(FLERR,"Illegal dump image command"); - bodyflag1 = force->numeric(FLERR,arg[iarg+2]); - bodyflag2 = force->numeric(FLERR,arg[iarg+3]); - iarg += 4; + } else if (strcmp(arg[iarg],"adiam") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal dump image command"); + adiam = NUMERIC; + adiamvalue = force->numeric(FLERR,arg[iarg+1]); + if (adiamvalue <= 0.0) error->all(FLERR,"Illegal dump image command"); + iarg += 2; } else if (strcmp(arg[iarg],"bond") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); @@ -174,6 +167,32 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : else error->all(FLERR,"Illegal dump image command"); iarg += 3; + + } else if (strcmp(arg[iarg],"line") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); + lineflag = YES; + if (strcmp(arg[iarg+1],"type") == 0) lcolor = TYPE; + else error->all(FLERR,"Illegal dump image command"); + ldiam = NUMERIC; + ldiamvalue = force->numeric(FLERR,arg[iarg+2]); + iarg += 3; + + } else if (strcmp(arg[iarg],"tri") == 0) { + if (iarg+2 > 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; + + } else if (strcmp(arg[iarg],"body") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal dump image command"); + bodyflag = YES; + if (strcmp(arg[iarg+1],"type") == 0) bodycolor = TYPE; + else error->all(FLERR,"Illegal dump image command"); + bodyflag1 = force->numeric(FLERR,arg[iarg+2]); + bodyflag2 = force->numeric(FLERR,arg[iarg+3]); + iarg += 4; + } else if (strcmp(arg[iarg],"size") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal dump image command"); int width = force->inumeric(FLERR,arg[iarg+1]); @@ -333,13 +352,26 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR,"Illegal dump image command"); } - // error check for bodyflag + // error checks and setup for lineflag, triflag, bodyflag - if (bodyflag) { - AtomVecBody *avec = (AtomVecBody *) atom->style_match("body"); - if (!avec) error->all(FLERR,"Dump image body yes requires atom style body"); - bptr = avec->bptr; + if (lineflag) { + avec_line = (AtomVecLine *) atom->style_match("line"); + if (!avec_line) + error->all(FLERR,"Dump image line requires atom style line"); } + if (triflag) { + avec_tri = (AtomVecTri *) atom->style_match("tri"); + if (!avec_tri) + error->all(FLERR,"Dump image tri requires atom style tri"); + } + if (bodyflag) { + avec_body = (AtomVecBody *) atom->style_match("body"); + if (!avec_body) + error->all(FLERR,"Dump image body yes requires atom style body"); + } + + extraflag = 0; + if (lineflag || triflag || bodyflag) extraflag = 1; // allocate image buffer now that image size is known @@ -679,18 +711,21 @@ void DumpImage::view_params() void DumpImage::create_image() { - int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus; + int i,j,k,m,n,itype,atom1,atom2,imol,iatom,btype,ibonus,drawflag; tagint tagprev; double diameter,delx,dely,delz; int *bodyvec; double **bodyarray; double *color,*color1,*color2; - double xmid[3]; + double xmid[3],pt1[3],pt2[3],pt3[3]; + double mat[3][3]; // render my atoms if (atomflag) { double **x = atom->x; + int *line = atom->line; + int *tri = atom->tri; int *body = atom->body; m = 0; @@ -719,18 +754,109 @@ void DumpImage::create_image() diameter = buf[m+1]; } - if (!body || !bodyflag || body[j] < 0) - image->draw_sphere(x[j],color,diameter); - else { - ibonus = body[i]; - n = bptr->image(ibonus,bodyflag1,bodyflag2,bodyvec,bodyarray); - for (k = 0; k < n; k++) { - if (bodyvec[k] == SPHERE) - image->draw_sphere(bodyarray[k],color,bodyarray[k][3]); - else if (bodyvec[k] == LINE) - image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3], - color,bodyarray[k][6],3); - } + // do not draw if line,tri,body keywords enabled and atom is one of those + + drawflag = 1; + if (extraflag) { + if (lineflag && line[j] >= 0) drawflag = 0; + if (triflag && tri[j] >= 0) drawflag = 0; + if (bodyflag && body[j] >= 0) drawflag = 0; + } + + if (drawflag) image->draw_sphere(x[j],color,diameter); + + m += size_one; + } + } + + // render atoms that are lines + + if (lineflag) { + double length,theta,dx,dy; + double **x = atom->x; + int *line = atom->line; + int *type = atom->type; + + for (i = 0; i < nchoose; i++) { + j = clist[i]; + if (line[j] < 0) continue; + + if (lcolor == TYPE) { + color = colortype[type[j]]; + } + + if (ldiam == NUMERIC) { + diameter = ldiamvalue; + } + + length = avec_line->bonus[line[j]].length; + theta = avec_line->bonus[line[j]].theta; + dx = 0.5*length*cos(theta); + dy = 0.5*length*sin(theta); + + pt1[0] = x[j][0] + dx; + pt1[1] = x[j][1] + dy; + pt1[2] = 0.0; + pt2[0] = x[j][0] - dx; + pt2[1] = x[j][1] - dy; + pt2[2] = 0.0; + + image->draw_cylinder(pt1,pt2,color,ldiamvalue,3); + } + } + + // render atoms that are triangles + + if (triflag) { + 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; + + if (tcolor == TYPE) { + color = colortype[type[j]]; + } + + MathExtra::quat_to_mat(avec_tri->bonus[tri[i]].quat,mat); + MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c1,pt1); + MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c2,pt2); + MathExtra::matvec(mat,avec_tri->bonus[tri[i]].c3,pt3); + MathExtra::add3(pt1,x[i],pt1); + MathExtra::add3(pt2,x[i],pt2); + MathExtra::add3(pt3,x[i],pt3); + + image->draw_triangle(pt1,pt2,pt3,color); + } + } + + // render atoms that are bodies + + if (bodyflag) { + Body *bptr = avec_body->bptr; + double **x = atom->x; + int *body = atom->body; + + m = 0; + for (i = 0; i < nchoose; i++) { + j = clist[i]; + if (body[j] < 0) continue; + + if (bodycolor == TYPE) { + itype = static_cast (buf[m]); + color = colortype[itype]; + } + + ibonus = body[i]; + n = bptr->image(ibonus,bodyflag1,bodyflag2,bodyvec,bodyarray); + for (k = 0; k < n; k++) { + if (bodyvec[k] == SPHERE) + image->draw_sphere(bodyarray[k],color,bodyarray[k][3]); + else if (bodyvec[k] == LINE) + image->draw_cylinder(&bodyarray[k][0],&bodyarray[k][3], + color,bodyarray[k][6],3); } m += size_one; diff --git a/src/dump_image.h b/src/dump_image.h index 9f3621d2c3..af6652c3ea 100644 --- a/src/dump_image.h +++ b/src/dump_image.h @@ -37,13 +37,24 @@ class DumpImage : public DumpCustom { int filetype; enum{PPM,JPG,PNG}; + int atomflag; // 0/1 for draw atoms int acolor,adiam; // what determines color/diam of atoms double adiamvalue; // atom diameter value - int atomflag,bondflag; // 0/1 for draw atoms,bonds + + int lineflag; // 0/1 for draw atoms as lines + 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 bodyflag; // 0/1 for draw atoms as bodies - double bodyflag1,bodyflag2; // user params for drawing bodies + int bodycolor; // what determines color of bodies + double bodyflag1,bodyflag2; // user-specified params for drawing bodies + + int bondflag; // 0/1 for draw bonds int bcolor,bdiam; // what determines color/diam of bonds double bdiamvalue; // bond diameter value + + int extraflag; // 0/1 for any of line/tri/body flag set char *thetastr,*phistr; // variables for view theta,phi int thetavar,phivar; // index to theta,phi vars int cflag; // static/dynamic box center @@ -64,8 +75,11 @@ class DumpImage : public DumpCustom { double *diamtype,*diamelement,*bdiamtype; // per-type diameters double **colortype,**colorelement,**bcolortype; // per-type colors + class AtomVecLine *avec_line; // ptrs to atom style (sub)classes + class AtomVecTri *avec_tri; + class AtomVecBody *avec_body; + class Image *image; // class that renders each image - class Body *bptr; // class for Body particles int *chooseghost; // extended choose array for comm double **bufcopy; // buffer for communicating bond/atom info diff --git a/src/fix_move.cpp b/src/fix_move.cpp index e1f522177b..dc7c21f3d2 100644 --- a/src/fix_move.cpp +++ b/src/fix_move.cpp @@ -26,7 +26,12 @@ #include "respa.h" #include "input.h" #include "variable.h" +#include "atom_vec_ellipsoid.h" +#include "atom_vec_line.h" +#include "atom_vec_tri.h" +#include "atom_vec_body.h" #include "math_const.h" +#include "math_extra.h" #include "memory.h" #include "error.h" @@ -37,6 +42,8 @@ using namespace MathConst; enum{LINEAR,WIGGLE,ROTATE,VARIABLE}; enum{EQUAL,ATOM}; +#define INERTIA 0.2 // moment of inertia prefactor for ellipsoid + /* ---------------------------------------------------------------------- */ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : @@ -216,7 +223,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : // set omega_rotate from period - if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = 2.0*MY_PI / period; + if (mstyle == WIGGLE || mstyle == ROTATE) omega_rotate = MY_2PI / period; // runit = unit vector along rotation axis @@ -229,24 +236,54 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : runit[2] = axis[2]/len; } - // set omega_flag if particles store omega - + // set flags for extra attributes particles may store + // relevant extra attributes = omega, angmom, theta, quat + omega_flag = atom->omega_flag; + angmom_flag = atom->angmom_flag; + radius_flag = atom->radius_flag; + ellipsoid_flag = atom->ellipsoid_flag; + line_flag = atom->line_flag; + tri_flag = atom->tri_flag; + body_flag = atom->body_flag; + + theta_flag = quat_flag = 0; + if (line_flag) theta_flag = 1; + if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1; + + extra_flag = 0; + if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1; + // perform initial allocation of atom-based array // register with Atom class xoriginal = NULL; + toriginal = NULL; + qoriginal = NULL; grow_arrays(atom->nmax); atom->add_callback(0); atom->add_callback(1); displace = velocity = NULL; + // AtomVec pointers to retrieve per-atom storage of extra quantities + + avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); + avec_line = (AtomVecLine *) atom->style_match("line"); + avec_tri = (AtomVecTri *) atom->style_match("tri"); + avec_body = (AtomVecBody *) atom->style_match("body"); + // xoriginal = initial unwrapped positions of atoms + // toriginal = initial theta of lines + // qoriginal = initial quat of extended particles double **x = atom->x; imageint *image = atom->image; + int *ellipsoid = atom->ellipsoid; + int *line = atom->line; + int *tri = atom->tri; + int *body = atom->body; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -255,6 +292,45 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) : else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0; } + if (theta_flag) { + for (int i = 0; i < nlocal; i++) { + if ((mask[i] & groupbit) && line[i] >= 0) + toriginal[i] = avec_line->bonus[line[i]].theta; + else toriginal[i] = 0.0; + } + } + + if (quat_flag) { + double *quat; + for (int i = 0; i < nlocal; i++) { + quat = NULL; + if (mask[i] & groupbit) { + if (ellipsoid_flag && ellipsoid[i] >= 0) + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (tri_flag && tri[i] >= 0) + quat = avec_tri->bonus[tri[i]].quat; + else if (body_flag && body[i] >= 0) + quat = avec_body->bonus[body[i]].quat; + } + if (quat) { + qoriginal[i][0] = quat[0]; + qoriginal[i][1] = quat[1]; + qoriginal[i][2] = quat[2]; + qoriginal[i][3] = quat[3]; + } else qoriginal[i][0] = qoriginal[i][1] = + qoriginal[i][2] = qoriginal[i][3] = 0.0; + } + } + + // nrestart = size of per-atom restart data + // nrestart = 1 + xorig + torig + qorig + + nrestart = 4; + if (theta_flag) nrestart++; + if (quat_flag) nrestart += 4; + + // time origin for movement = current timestep + time_origin = update->ntimestep; } @@ -270,6 +346,8 @@ FixMove::~FixMove() // delete locally stored arrays memory->destroy(xoriginal); + memory->destroy(toriginal); + memory->destroy(qoriginal); memory->destroy(displace); memory->destroy(velocity); @@ -381,9 +459,12 @@ void FixMove::init() void FixMove::initial_integrate(int vflag) { - double dtfm; - double xold[3],a[3],b[3],c[3],d[3],disp[3]; + int flag; double ddotr,dx,dy,dz; + double dtfm,theta_new; + double xold[3],a[3],b[3],c[3],d[3],disp[3],w[3],ex[3],ey[3],ez[3]; + double inertia_ellipsoid[3],qrotate[4]; + double *quat,*inertia,*shape; double delta = (update->ntimestep - time_origin) * dt; @@ -391,10 +472,17 @@ void FixMove::initial_integrate(int vflag) double **v = atom->v; double **f = atom->f; double **omega = atom->omega; + double **angmom = atom->angmom; + double *radius = atom->radius; double *rmass = atom->rmass; double *mass = atom->mass; int *type = atom->type; + int *ellipsoid = atom->ellipsoid; + int *line = atom->line; + int *tri = atom->tri; + int *body = atom->body; int *mask = atom->mask; + int nlocal = atom->nlocal; // for linear: X = X0 + V*dt @@ -553,17 +641,89 @@ void FixMove::initial_integrate(int vflag) v[i][0] = omega_rotate * (runit[1]*disp[2] - runit[2]*disp[1]); v[i][1] = omega_rotate * (runit[2]*disp[0] - runit[0]*disp[2]); v[i][2] = omega_rotate * (runit[0]*disp[1] - runit[1]*disp[0]); - if (omega_flag) { - omega[i][0] = omega_rotate*runit[0]; - omega[i][1] = omega_rotate*runit[1]; - omega[i][2] = omega_rotate*runit[2]; - } + + // set any extra attributes affected by rotation + + if (extra_flag) { + + // omega for spheres and lines + + if (omega_flag) { + flag = 0; + if (radius_flag && radius[i] > 0.0) flag = 1; + if (line_flag && line[i] >= 0.0) flag = 1; + if (flag) { + omega[i][0] = omega_rotate*runit[0]; + omega[i][1] = omega_rotate*runit[1]; + omega[i][2] = omega_rotate*runit[2]; + } + } + + // angmom for ellipsoids, tris, and bodies + + if (angmom_flag) { + quat = inertia = NULL; + if (ellipsoid_flag && ellipsoid[i] >= 0) { + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + shape = avec_ellipsoid->bonus[ellipsoid[i]].shape; + inertia_ellipsoid[0] = + INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]); + inertia_ellipsoid[1] = + INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]); + inertia_ellipsoid[2] = + INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]); + inertia = inertia_ellipsoid; + } else if (tri_flag && tri[i] >= 0) { + quat = avec_tri->bonus[tri[i]].quat; + inertia = avec_tri->bonus[tri[i]].inertia; + } else if (body_flag && body[i] >= 0) { + quat = avec_body->bonus[body[i]].quat; + inertia = avec_body->bonus[body[i]].inertia; + } + if (quat) { + w[0] = omega_rotate*runit[0]; + w[1] = omega_rotate*runit[1]; + w[2] = omega_rotate*runit[2]; + MathExtra::q_to_exyz(quat,ex,ey,ez); + MathExtra::omega_to_angmom(w,ex,ey,ez,inertia,angmom[i]); + } + } + + // theta for lines + + if (theta_flag && line[i] >= 0.0) { + theta_new = fmod(toriginal[i]+arg,MY_2PI); + avec_line->bonus[atom->line[i]].theta = theta_new; + } + + // quats for ellipsoids, tris, and bodies + + if (quat_flag) { + quat = NULL; + if (ellipsoid_flag && ellipsoid[i] >= 0) + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (tri_flag && tri[i] >= 0) + 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); + } + } + } domain->remap_near(x[i],xold); } } // for variable: compute x,v from variables + // NOTE: also allow for changes to extra attributes? + // omega, angmom, theta, quat + // only necessary if prescribed motion involves rotation } else if (mstyle == VARIABLE) { @@ -809,6 +969,8 @@ void FixMove::final_integrate_respa(int ilevel, int iloop) double FixMove::memory_usage() { double bytes = atom->nmax*3 * sizeof(double); + if (theta_flag) bytes += atom->nmax * sizeof(double); + if (quat_flag) bytes += atom->nmax*4 * sizeof(double); if (displaceflag) bytes += atom->nmax*3 * sizeof(double); if (velocityflag) bytes += atom->nmax*3 * sizeof(double); return bytes; @@ -850,6 +1012,8 @@ void FixMove::restart(char *buf) void FixMove::grow_arrays(int nmax) { memory->grow(xoriginal,nmax,3,"move:xoriginal"); + if (theta_flag) memory->grow(toriginal,nmax,"move:toriginal"); + if (quat_flag) memory->grow(qoriginal,nmax,4,"move:qoriginal"); array_atom = xoriginal; } @@ -862,6 +1026,13 @@ void FixMove::copy_arrays(int i, int j, int delflag) xoriginal[j][0] = xoriginal[i][0]; xoriginal[j][1] = xoriginal[i][1]; xoriginal[j][2] = xoriginal[i][2]; + if (theta_flag) toriginal[j] = toriginal[i]; + if (quat_flag) { + qoriginal[j][0] = qoriginal[i][0]; + qoriginal[j][1] = qoriginal[i][1]; + qoriginal[j][2] = qoriginal[i][2]; + qoriginal[j][3] = qoriginal[i][3]; + } } /* ---------------------------------------------------------------------- @@ -870,8 +1041,15 @@ void FixMove::copy_arrays(int i, int j, int delflag) void FixMove::set_arrays(int i) { + double theta; + double *quat; + double **x = atom->x; imageint *image = atom->image; + int *ellipsoid = atom->ellipsoid; + int *line = atom->line; + int *tri = atom->tri; + int *body = atom->body; int *mask = atom->mask; // particle not in group @@ -932,6 +1110,33 @@ void FixMove::set_arrays(int i) xoriginal[i][0] = point[0] + c[0] + disp[0]; xoriginal[i][1] = point[1] + c[1] + disp[1]; xoriginal[i][2] = point[2] + c[2] + disp[2]; + + // set theta and quat extra attributes affected by rotation + + if (extra_flag) { + + // theta for lines + + if (theta_flag && line[i] >= 0.0) { + theta = avec_line->bonus[atom->line[i]].theta; + toriginal[i] = theta - 0.0; // NOTE: edit this line + } + + // quats for ellipsoids, tris, and bodies + + if (quat_flag) { + quat = NULL; + if (ellipsoid_flag && ellipsoid[i] >= 0) + quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; + else if (tri_flag && tri[i] >= 0) + quat = avec_tri->bonus[tri[i]].quat; + else if (body_flag && body[i] >= 0) + quat = avec_body->bonus[body[i]].quat; + if (quat) { + // qoriginal = f(quat,-delta); // NOTE: edit this line + } + } + } } } @@ -941,10 +1146,18 @@ void FixMove::set_arrays(int i) int FixMove::pack_exchange(int i, double *buf) { - buf[0] = xoriginal[i][0]; - buf[1] = xoriginal[i][1]; - buf[2] = xoriginal[i][2]; - return 3; + int n = 0; + buf[n++] = xoriginal[i][0]; + buf[n++] = xoriginal[i][1]; + buf[n++] = xoriginal[i][2]; + if (theta_flag) buf[n++] = toriginal[i]; + if (quat_flag) { + buf[n++] = qoriginal[i][0]; + buf[n++] = qoriginal[i][1]; + buf[n++] = qoriginal[i][2]; + buf[n++] = qoriginal[i][3]; + } + return n; } /* ---------------------------------------------------------------------- @@ -953,10 +1166,18 @@ int FixMove::pack_exchange(int i, double *buf) int FixMove::unpack_exchange(int nlocal, double *buf) { - xoriginal[nlocal][0] = buf[0]; - xoriginal[nlocal][1] = buf[1]; - xoriginal[nlocal][2] = buf[2]; - return 3; + int n = 0; + xoriginal[nlocal][0] = buf[n++]; + xoriginal[nlocal][1] = buf[n++]; + xoriginal[nlocal][2] = buf[n++]; + if (theta_flag) toriginal[nlocal] = buf[n++]; + if (quat_flag) { + qoriginal[nlocal][0] = buf[n++]; + qoriginal[nlocal][1] = buf[n++]; + qoriginal[nlocal][2] = buf[n++]; + qoriginal[nlocal][3] = buf[n++]; + } + return n; } /* ---------------------------------------------------------------------- @@ -965,11 +1186,19 @@ int FixMove::unpack_exchange(int nlocal, double *buf) int FixMove::pack_restart(int i, double *buf) { - buf[0] = 4; - buf[1] = xoriginal[i][0]; - buf[2] = xoriginal[i][1]; - buf[3] = xoriginal[i][2]; - return 4; + int n = 1; + buf[n++] = xoriginal[i][0]; + buf[n++] = xoriginal[i][1]; + buf[n++] = xoriginal[i][2]; + if (theta_flag) buf[n++] = toriginal[i]; + if (quat_flag) { + buf[n++] = qoriginal[i][0]; + buf[n++] = qoriginal[i][1]; + buf[n++] = qoriginal[i][2]; + buf[n++] = qoriginal[i][3]; + } + buf[0] = n; + return n; } /* ---------------------------------------------------------------------- @@ -989,6 +1218,13 @@ void FixMove::unpack_restart(int nlocal, int nth) xoriginal[nlocal][0] = extra[nlocal][m++]; xoriginal[nlocal][1] = extra[nlocal][m++]; xoriginal[nlocal][2] = extra[nlocal][m++]; + if (theta_flag) toriginal[nlocal] = extra[nlocal][m++]; + if (quat_flag) { + qoriginal[nlocal][0] = extra[nlocal][m++]; + qoriginal[nlocal][1] = extra[nlocal][m++]; + qoriginal[nlocal][2] = extra[nlocal][m++]; + qoriginal[nlocal][3] = extra[nlocal][m++]; + } } /* ---------------------------------------------------------------------- @@ -997,7 +1233,7 @@ void FixMove::unpack_restart(int nlocal, int nth) int FixMove::maxsize_restart() { - return 4; + return nrestart; } /* ---------------------------------------------------------------------- @@ -1006,7 +1242,7 @@ int FixMove::maxsize_restart() int FixMove::size_restart(int nlocal) { - return 4; + return nrestart; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_move.h b/src/fix_move.h index bcd8775708..977ec88eae 100644 --- a/src/fix_move.h +++ b/src/fix_move.h @@ -61,13 +61,23 @@ class FixMove : public Fix { double dt,dtv,dtf; int xvar,yvar,zvar,vxvar,vyvar,vzvar; int xvarstyle,yvarstyle,zvarstyle,vxvarstyle,vyvarstyle,vzvarstyle; - int omega_flag,nlevels_respa; + int extra_flag,omega_flag,angmom_flag; + int radius_flag,ellipsoid_flag,line_flag,tri_flag,body_flag; + int theta_flag,quat_flag; + int nlevels_respa,nrestart; int time_origin; double **xoriginal; // original coords of atoms + double *toriginal; // original theta of atoms + double **qoriginal; // original quat of atoms int displaceflag,velocityflag; int maxatom; double **displace,**velocity; + + class AtomVecEllipsoid *avec_ellipsoid; + class AtomVecLine *avec_line; + class AtomVecTri *avec_tri; + class AtomVecBody *avec_body; }; } diff --git a/src/math_extra.cpp b/src/math_extra.cpp index 0b33fc59c1..7095412b3e 100644 --- a/src/math_extra.cpp +++ b/src/math_extra.cpp @@ -230,6 +230,7 @@ void richardson(double *q, double *m, double *w, double *moments, double dtq) apply evolution operators to quat, quat momentum Miller et al., J Chem Phys. 116, 8649-8659 (2002) ------------------------------------------------------------------------- */ + void no_squish_rotate(int k, double *p, double *q, double *inertia, double dt) { diff --git a/src/math_extra.h b/src/math_extra.h index f768715ca0..58cae1df28 100755 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -76,7 +76,8 @@ namespace MathExtra { void rotate(double matrix[3][3], int i, int j, int k, int l, double s, double tau); void richardson(double *q, double *m, double *w, double *moments, double dtq); - void no_squish_rotate(int k, double *p, double *q, double *inertia, double dt); + void no_squish_rotate(int k, double *p, double *q, double *inertia, + double dt); // shape matrix operations // upper-triangular 3x3 matrix stored in Voigt notation as 6-vector @@ -151,7 +152,8 @@ inline void MathExtra::normalize3(const double *v, double *ans) scale a vector to length ------------------------------------------------------------------------- */ -inline void MathExtra::snormalize3(const double length, const double *v, double *ans) +inline void MathExtra::snormalize3(const double length, const double *v, + double *ans) { double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); ans[0] = v[0]*scale; @@ -312,7 +314,7 @@ void MathExtra::times3_diag(const double m[3][3], const double *d, ------------------------------------------------------------------------- */ inline void MathExtra::plus3(const double m[3][3], const double m2[3][3], - double ans[3][3]) + double ans[3][3]) { ans[0][0] = m[0][0]+m2[0][0]; ans[0][1] = m[0][1]+m2[0][1]; @@ -330,7 +332,7 @@ inline void MathExtra::plus3(const double m[3][3], const double m2[3][3], ------------------------------------------------------------------------- */ inline void MathExtra::times3(const double m[3][3], const double m2[3][3], - double ans[3][3]) + double ans[3][3]) { ans[0][0] = m[0][0]*m2[0][0] + m[0][1]*m2[1][0] + m[0][2]*m2[2][0]; ans[0][1] = m[0][0]*m2[0][1] + m[0][1]*m2[1][1] + m[0][2]*m2[2][1]; @@ -601,7 +603,8 @@ inline void MathExtra::axisangle_to_quat(const double *v, const double angle, Apply principal rotation generator about x to rotation matrix m ------------------------------------------------------------------------- */ -inline void MathExtra::rotation_generator_x(const double m[3][3], double ans[3][3]) +inline void MathExtra::rotation_generator_x(const double m[3][3], + double ans[3][3]) { ans[0][0] = 0; ans[0][1] = -m[0][2];