From ca8dfd3d7c07ac0560cfe8230e705cc068fa1716 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 20 Oct 2011 14:56:49 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7146
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/SRD/fix_srd.cpp | 1162 +++++++++++++++++++++++----------
src/SRD/fix_srd.h | 52 +-
src/atom.cpp | 12 +-
src/atom.h | 4 +-
src/compute_property_atom.cpp | 428 +++++++++++-
src/compute_property_atom.h | 19 +-
src/fix_rigid.cpp | 325 ++++++---
src/fix_rigid.h | 11 +-
src/read_data.cpp | 99 ++-
src/read_data.h | 6 +-
src/set.cpp | 136 +++-
11 files changed, 1687 insertions(+), 567 deletions(-)
diff --git a/src/SRD/fix_srd.cpp b/src/SRD/fix_srd.cpp
index 411beb9c9d..b673b9b7b5 100644
--- a/src/SRD/fix_srd.cpp
+++ b/src/SRD/fix_srd.cpp
@@ -22,6 +22,8 @@
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
#include "group.h"
#include "update.h"
#include "force.h"
@@ -42,7 +44,7 @@ using namespace LAMMPS_NS;
using namespace MathConst;
enum{SLIP,NOSLIP};
-enum{SPHERE,ELLIPSOID,WALL};
+enum{SPHERE,ELLIPSOID,LINE,TRIANGLE,WALL};
enum{INSIDE_ERROR,INSIDE_WARN,INSIDE_IGNORE};
enum{BIG_MOVE,SRD_MOVE,SRD_ROTATE};
enum{CUBIC_ERROR,CUBIC_WARN};
@@ -50,10 +52,13 @@ enum{SHIFT_NO,SHIFT_YES,SHIFT_POSSIBLE};
enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp
-#define INERTIA 0.4
-#define ATOMPERBIN 10
+#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
+
+#define ATOMPERBIN 30
#define BIG 1.0e20
#define VBINSIZE 5
+#define TOLERANCE 0.00001
+#define MAXITER 20
//#define SRD_DEBUG 1
//#define SRD_DEBUG_ATOMID 58
@@ -95,7 +100,7 @@ FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
cubictol = 0.01;
shiftuser = SHIFT_NO;
shiftseed = 0;
- streamflag = 1;
+ tstat = 0;
int iarg = 8;
while (iarg < narg) {
@@ -156,10 +161,10 @@ FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
else error->all(FLERR,"Illegal fix srd command");
shiftseed = atoi(arg[iarg+2]);
iarg += 3;
- } else if (strcmp(arg[iarg],"stream") == 0) {
+ } else if (strcmp(arg[iarg],"tstat") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix srd command");
- if (strcmp(arg[iarg+1],"yes") == 0) streamflag = 1;
- else if (strcmp(arg[iarg+1],"no") == 0) streamflag = 0;
+ if (strcmp(arg[iarg+1],"no") == 0) tstat = 0;
+ else if (strcmp(arg[iarg+1],"yes") == 0) tstat = 1;
else error->all(FLERR,"Illegal fix srd command");
iarg += 2;
} else error->all(FLERR,"Illegal fix srd command");
@@ -168,7 +173,8 @@ FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// error check
if (nevery <= 0) error->all(FLERR,"Illegal fix srd command");
- if (bigexist && biggroup < 0) error->all(FLERR,"Could not find fix srd group ID");
+ if (bigexist && biggroup < 0)
+ error->all(FLERR,"Could not find fix srd group ID");
if (gridsrd <= 0.0) error->all(FLERR,"Illegal fix srd command");
if (temperature_srd <= 0.0) error->all(FLERR,"Illegal fix srd command");
if (seed <= 0) error->all(FLERR,"Illegal fix srd command");
@@ -176,7 +182,8 @@ FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
if (maxbounceallow < 0) error->all(FLERR,"Illegal fix srd command");
if (lamdaflag && lamda <= 0.0) error->all(FLERR,"Illegal fix srd command");
if (gridsearch <= 0.0) error->all(FLERR,"Illegal fix srd command");
- if (cubictol < 0.0 || cubictol > 1.0) error->all(FLERR,"Illegal fix srd command");
+ if (cubictol < 0.0 || cubictol > 1.0)
+ error->all(FLERR,"Illegal fix srd command");
if ((shiftuser == SHIFT_YES || shiftuser == SHIFT_POSSIBLE) &&
shiftseed <= 0) error->all(FLERR,"Illegal fix srd command");
@@ -236,6 +243,8 @@ FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// atom style pointers to particles that store bonus info
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
// fix parameters
@@ -289,7 +298,8 @@ void FixSRD::init()
{
// error checks
- if (force->newton_pair == 0) error->all(FLERR,"Fix srd requires newton pair on");
+ if (force->newton_pair == 0)
+ error->all(FLERR,"Fix srd requires newton pair on");
if (bigexist && comm->ghost_velocity == 0)
error->all(FLERR,"Fix srd requires ghost atoms store velocity");
if (bigexist && collidestyle == NOSLIP && !atom->torque_flag)
@@ -319,19 +329,21 @@ void FixSRD::init()
fwall = wallfix->fwall;
walltrigger = 0.5 * neighbor->skin;
if (wallfix->overlap && overlap == 0 && me == 0)
- error->warning(FLERR,"Fix SRD walls overlap but fix srd overlap not set");
+ error->warning(FLERR,
+ "Fix SRD walls overlap but fix srd overlap not set");
}
}
// set change_flags if box size or shape changes
- change_size = change_shape = 0;
+ change_size = change_shape = deformflag = 0;
if (domain->nonperiodic == 2) change_size = 1;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->box_change) {
if (modify->fix[i]->box_change_size) change_size = 1;
if (modify->fix[i]->box_change_shape) change_shape = 1;
if (strcmp(modify->fix[i]->style,"deform") == 0) {
+ deformflag = 1;
FixDeform *deform = (FixDeform *) modify->fix[i];
if (deform->box_change_shape && deform->remapflag != V_REMAP)
error->all(FLERR,"Using fix srd with inconsistent "
@@ -339,6 +351,10 @@ void FixSRD::init()
}
}
+ if (deformflag && tstat == 0 && me == 0)
+ error->warning(FLERR,
+ "Using fix srd with box deformation but no SRD thermostat");
+
// parameterize based on current box volume
dimension = domain->dimension;
@@ -391,7 +407,8 @@ void FixSRD::setup(int vflag)
setup_bounds();
if (dist_srd_reneigh < nevery*dt_big*vmax && me == 0)
- error->warning(FLERR,"Fix srd SRD moves may trigger frequent reneighboring");
+ error->warning(FLERR,
+ "Fix srd SRD moves may trigger frequent reneighboring");
// setup search bins and search stencil based on these distances
@@ -407,6 +424,7 @@ void FixSRD::setup(int vflag)
reneighflag = BIG_MOVE;
pre_neighbor();
}
+
/* ----------------------------------------------------------------------
assign SRD particles to bins
assign big particles to all bins they overlap
@@ -416,6 +434,7 @@ void FixSRD::pre_neighbor()
{
int i,j,m,ix,iy,iz,jx,jy,jz,ibin,jbin,lo,hi;
double rsq,cutbinsq;
+ double xlamda[3];
// grow SRD per-atom bin arrays if necessary
@@ -561,7 +580,8 @@ void FixSRD::pre_neighbor()
if (side == 0) {
hi = static_cast<int> ((xwall[m]+delta-xblo2)*bininv2x);
if (hi < 0) continue;
- if (hi >= nbin2x) error->all(FLERR,"Fix SRD: bad search bin assignment");
+ if (hi >= nbin2x) error->all(FLERR,
+ "Fix SRD: bad search bin assignment");
lo = 0;
} else {
lo = static_cast<int> ((xwall[m]-delta-xblo2)*bininv2x);
@@ -583,7 +603,8 @@ void FixSRD::pre_neighbor()
if (side == 0) {
hi = static_cast<int> ((xwall[m]+delta-yblo2)*bininv2y);
if (hi < 0) continue;
- if (hi >= nbin2y) error->all(FLERR,"Fix SRD: bad search bin assignment");
+ if (hi >= nbin2y) error->all(FLERR,
+ "Fix SRD: bad search bin assignment");
lo = 0;
} else {
lo = static_cast<int> ((xwall[m]-delta-yblo2)*bininv2y);
@@ -605,7 +626,8 @@ void FixSRD::pre_neighbor()
if (side == 0) {
hi = static_cast<int> ((xwall[m]+delta-zblo2)*bininv2z);
if (hi < 0) continue;
- if (hi >= nbin2z) error->all(FLERR,"Fix SRD: bad search bin assignment");
+ if (hi >= nbin2z) error->all(FLERR,
+ "Fix SRD: bad search bin assignment");
lo = 0;
} else {
lo = static_cast<int> ((xwall[m]-delta-zblo2)*bininv2z);
@@ -645,6 +667,7 @@ void FixSRD::pre_neighbor()
void FixSRD::post_force(int vflag)
{
int i,m,ix,iy,iz;
+ double xlamda[3];
// zero per-timestep stats
@@ -765,19 +788,20 @@ void FixSRD::post_force(int vflag)
/* ----------------------------------------------------------------------
reset SRD velocities
- may perform random shifting by 1/2 bin in each dimension
+ may perform random shifting by up to 1/2 bin in each dimension
called at pre-neighbor stage when all SRDs are now inside my sub-domain
- for triclinic, will set mean velocity to box deformation velocity
+ if tstat, then thermostat SRD particles as well, including streaming effects
------------------------------------------------------------------------- */
void FixSRD::reset_velocities()
{
int i,j,n,ix,iy,iz,ibin,axis,sign,irandom;
- double u[3],vave[3];
- double vx,vy,vz,vsq;
- double *vold,*vnew,*xlamda;
+ double u[3],vsum[3];
+ double vx,vy,vz,vsq,tbin,scale;
+ double *vave,*vnew,*xlamda;
+ double vstream[3];
- // if requested, perform a dynamic shift
+ // if requested, perform a dynamic shift of bin positions
if (shiftflag) {
double *boxlo;
@@ -831,39 +855,23 @@ void FixSRD::reset_velocities()
if (triclinic) domain->lamda2x(nlocal);
- if (triclinic && streamflag) {
- double *h_rate = domain->h_rate;
- double *h_ratelo = domain->h_ratelo;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
- xlamda = x[i];
- v[i][0] -= h_rate[0]*xlamda[0] + h_rate[5]*xlamda[1] +
- h_rate[4]*xlamda[2] + h_ratelo[0];
- v[i][1] -= h_rate[1]*xlamda[1] + h_rate[3]*xlamda[2] + h_ratelo[1];
- v[i][2] -= h_rate[2]*xlamda[2] + h_ratelo[2];
- }
- }
-
// for each bin I have particles contributing to:
- // compute ave v and v^2 of particles in that bin
+ // compute summed v of particles in that bin
// if I own the bin, set its random value, else set to 0.0
for (i = 0; i < nbins; i++) {
n = 0;
- vave[0] = vave[1] = vave[2] = 0.0;
+ vsum[0] = vsum[1] = vsum[2] = 0.0;
for (j = binhead[i]; j >= 0; j = binnext[j]) {
- vx = v[j][0];
- vy = v[j][1];
- vz = v[j][2];
- vave[0] += vx;
- vave[1] += vy;
- vave[2] += vz;
+ vsum[0] += v[j][0];
+ vsum[1] += v[j][1];
+ vsum[2] += v[j][2];
n++;
}
- vbin[i].vave[0] = vave[0];
- vbin[i].vave[1] = vave[1];
- vbin[i].vave[2] = vave[2];
+ vbin[i].vsum[0] = vsum[0];
+ vbin[i].vsum[1] = vsum[1];
+ vbin[i].vsum[2] = vsum[2];
vbin[i].n = n;
if (vbin[i].owner) vbin[i].random = random->uniform();
else vbin[i].random = 0.0;
@@ -874,24 +882,43 @@ void FixSRD::reset_velocities()
if (shifts[shiftflag].commflag) vbin_comm(shiftflag);
// for each bin I have particles contributing to:
- // reassign particle velocity by rotation around a random axis
- // accumulate T_srd for each bin I own
- // for triclinic, replace mean velocity with stream velocity
+ // compute vave over particles in bin
+ // thermal velocity of each particle = v - vave
+ // rotate thermal vel of each particle around one of 6 random axes
+ // add vave back to each particle
+ // thermostat if requested:
+ // if no deformation, rescale thermal vel to temperature
+ // if deformation, rescale thermal vel and change vave to vstream
+ // these are settings for extra dof_temp, dof_tstat to subtract
+ // (not sure why these settings work best)
+ // no deformation, no tstat: dof_temp = 1
+ // yes deformation, no tstat: doesn't matter, system will not equilibrate
+ // no deformation, yes tstat: dof_temp = dof_tstat = 1
+ // yes deformation, yes tstat: dof_temp = dof_tstat = 0
+ // accumulate final T_srd for each bin I own
+
+ double tfactor = force->mvv2e * mass_srd / (dimension * force->boltz);
+ int dof_temp = 1;
+ int dof_tstat;
+ if (tstat) {
+ if (deformflag) dof_tstat = dof_temp = 0;
+ else dof_tstat = 1;
+ }
srd_bin_temp = 0.0;
srd_bin_count = 0;
if (dimension == 2) axis = 2;
+ double *h_rate = domain->h_rate;
+ double *h_ratelo = domain->h_ratelo;
for (i = 0; i < nbins; i++) {
n = vbin[i].n;
if (n == 0) continue;
- vold = vbin[i].vave;
- vold[0] /= n;
- vold[1] /= n;
- vold[2] /= n;
-
- vnew = vold;
+ vave = vbin[i].vsum;
+ vave[0] /= n;
+ vave[1] /= n;
+ vave[2] /= n;
irandom = static_cast<int> (6.0*vbin[i].random);
sign = irandom % 2;
@@ -900,48 +927,64 @@ void FixSRD::reset_velocities()
vsq = 0.0;
for (j = binhead[i]; j >= 0; j = binnext[j]) {
if (axis == 0) {
- u[0] = v[j][0]-vold[0];
- u[1] = sign ? v[j][2]-vold[2] : vold[2]-v[j][2];
- u[2] = sign ? vold[1]-v[j][1] : v[j][1]-vold[1];
+ u[0] = v[j][0]-vave[0];
+ u[1] = sign ? v[j][2]-vave[2] : vave[2]-v[j][2];
+ u[2] = sign ? vave[1]-v[j][1] : v[j][1]-vave[1];
} else if (axis == 1) {
- u[1] = v[j][1]-vold[1];
- u[0] = sign ? v[j][2]-vold[2] : vold[2]-v[j][2];
- u[2] = sign ? vold[0]-v[j][0] : v[j][0]-vold[0];
+ u[1] = v[j][1]-vave[1];
+ u[0] = sign ? v[j][2]-vave[2] : vave[2]-v[j][2];
+ u[2] = sign ? vave[0]-v[j][0] : v[j][0]-vave[0];
} else {
- u[2] = v[j][2]-vold[2];
- u[1] = sign ? v[j][0]-vold[0] : vold[0]-v[j][0];
- u[0] = sign ? vold[1]-v[j][1] : v[j][1]-vold[1];
+ u[2] = v[j][2]-vave[2];
+ u[1] = sign ? v[j][0]-vave[0] : vave[0]-v[j][0];
+ u[0] = sign ? vave[1]-v[j][1] : v[j][1]-vave[1];
}
vsq += u[0]*u[0] + u[1]*u[1] + u[2]*u[2];
- v[j][0] = u[0] + vnew[0];
- v[j][1] = u[1] + vnew[1];
- v[j][2] = u[2] + vnew[2];
+ v[j][0] = u[0] + vave[0];
+ v[j][1] = u[1] + vave[1];
+ v[j][2] = u[2] + vave[2];
}
- // sum partial contribution of my particles to T even if I don't own bin
- // but only count bin if I own it, so that bin is counted exactly once
+ if (tstat && n > 1) {
+ if (deformflag) {
+ xlamda = vbin[i].xctr;
+ vstream[0] = h_rate[0]*xlamda[0] + h_rate[5]*xlamda[1] +
+ h_rate[4]*xlamda[2] + h_ratelo[0];
+ vstream[1] = h_rate[1]*xlamda[1] + h_rate[3]*xlamda[2] + h_ratelo[1];
+ vstream[2] = h_rate[2]*xlamda[2] + h_ratelo[2];
+ } else {
+ vstream[0] = vave[0];
+ vstream[1] = vave[1];
+ vstream[2] = vave[2];
+ }
- if (n > 1) {
- srd_bin_temp += vsq / (n-1);
- if (vbin[i].owner) srd_bin_count++;
+ // tbin = thermal temperature of particles in bin
+ // scale = scale factor for thermal velocity
+
+ tbin = vsq/(n-dof_tstat) * tfactor;
+ scale = sqrt(temperature_srd/tbin);
+
+ vsq = 0.0;
+ for (j = binhead[i]; j >= 0; j = binnext[j]) {
+ u[0] = (v[j][0] - vave[0]) * scale;
+ u[1] = (v[j][1] - vave[1]) * scale;
+ u[2] = (v[j][2] - vave[2]) * scale;
+ vsq += u[0]*u[0] + u[1]*u[1] + u[2]*u[2];
+ v[j][0] = u[0] + vstream[0];
+ v[j][1] = u[1] + vstream[1];
+ v[j][2] = u[2] + vstream[2];
+ }
}
- }
- srd_bin_temp *= force->mvv2e * mass_srd / (dimension * force->boltz);
+ // sum partial contribution of my particles to T even if I don't own bin
+ // but only count bin if I own it, so each bin is counted exactly once
- if (triclinic && streamflag) {
- double *h_rate = domain->h_rate;
- double *h_ratelo = domain->h_ratelo;
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
- xlamda = x[i];
- v[i][0] += h_rate[0]*xlamda[0] + h_rate[5]*xlamda[1] +
- h_rate[4]*xlamda[2] + h_ratelo[0];
- v[i][1] += h_rate[1]*xlamda[1] + h_rate[3]*xlamda[2] + h_ratelo[1];
- v[i][2] += h_rate[2]*xlamda[2] + h_ratelo[2];
- }
+ if (n > 1) srd_bin_temp += vsq/(n-dof_temp);
+ if (vbin[i].owner) srd_bin_count++;
}
+ srd_bin_temp *= tfactor;
+
// rescale any too-large velocities
for (i = 0; i < nlocal; i++)
@@ -1032,9 +1075,9 @@ void FixSRD::vbin_pack(BinAve *vbin, int n, int *list, double *buf)
for (int i = 0; i < n; i++) {
j = list[i];
buf[m++] = vbin[j].n;
- buf[m++] = vbin[j].vave[0];
- buf[m++] = vbin[j].vave[1];
- buf[m++] = vbin[j].vave[2];
+ buf[m++] = vbin[j].vsum[0];
+ buf[m++] = vbin[j].vsum[1];
+ buf[m++] = vbin[j].vsum[2];
buf[m++] = vbin[j].random;
}
}
@@ -1050,9 +1093,9 @@ void FixSRD::vbin_unpack(double *buf, BinAve *vbin, int n, int *list)
for (int i = 0; i < n; i++) {
j = list[i];
vbin[j].n += static_cast<int> (buf[m++]);
- vbin[j].vave[0] += buf[m++];
- vbin[j].vave[1] += buf[m++];
- vbin[j].vave[2] += buf[m++];
+ vbin[j].vsum[0] += buf[m++];
+ vbin[j].vsum[1] += buf[m++];
+ vbin[j].vsum[2] += buf[m++];
vbin[j].random += buf[m++];
}
}
@@ -1065,7 +1108,7 @@ void FixSRD::vbin_unpack(double *buf, BinAve *vbin, int n, int *list)
void FixSRD::collisions_single()
{
- int i,j,k,m,type,mbig,ibin,ibounce,inside,collide_flag;
+ int i,j,k,m,type,nbig,ibin,ibounce,inside,collide_flag,lineside;
double dt,t_remain;
double norm[3],xscoll[3],xbcoll[3],vsnew[3];
Big *big;
@@ -1097,11 +1140,11 @@ void FixSRD::collisions_single()
dt = dt_big;
while (collide_flag) {
- mbig = nbinbig[ibin];
- if (ibounce == 0) ncheck += mbig;
+ nbig = nbinbig[ibin];
+ if (ibounce == 0) ncheck += nbig;
collide_flag = 0;
- for (m = 0; m < mbig; m++) {
+ for (m = 0; m < nbig; m++) {
k = binbig[ibin][m];
big = &biglist[k];
j = big->index;
@@ -1161,17 +1204,11 @@ void FixSRD::collisions_single()
}
if (collidestyle == SLIP) {
- if (type == SPHERE)
- slip_sphere(v[i],v[j],norm,vsnew);
- else if (type == ELLIPSOID)
- slip_ellipsoid(v[i],v[j],x[j],big,xscoll,norm,vsnew);
- else
- slip_wall(v[i],j,norm,vsnew);
+ if (type != WALL) slip(v[i],v[j],x[j],big,xscoll,norm,vsnew);
+ else slip_wall(v[i],j,norm,vsnew);
} else {
- if (type != WALL)
- noslip(v[i],v[j],x[j],big,xscoll,norm,vsnew);
- else
- noslip_wall(v[i],j,xscoll,norm,vsnew);
+ if (type != WALL) noslip(v[i],v[j],x[j],big,-1, xscoll,norm,vsnew);
+ else noslip(v[i],NULL,x[j],big,j,xscoll,norm,vsnew);
}
if (dimension == 2) vsnew[2] = 0.0;
@@ -1222,7 +1259,7 @@ void FixSRD::collisions_single()
void FixSRD::collisions_multi()
{
- int i,j,k,m,type,mbig,ibin,ibounce,inside,jfirst,typefirst;
+ int i,j,k,m,type,nbig,ibin,ibounce,inside,jfirst,typefirst,jlast;
double dt,t_remain,t_first;
double norm[3],xscoll[3],xbcoll[3],vsnew[3];
double normfirst[3],xscollfirst[3],xbcollfirst[3];
@@ -1251,22 +1288,31 @@ void FixSRD::collisions_multi()
if (nbinbig[ibin] == 0) continue;
ibounce = 0;
+ jlast = -1;
dt = dt_big;
while (1) {
- mbig = nbinbig[ibin];
- if (ibounce == 0) ncheck += mbig;
+ nbig = nbinbig[ibin];
+ if (ibounce == 0) ncheck += nbig;
t_first = 0.0;
- for (m = 0; m < mbig; m++) {
+ for (m = 0; m < nbig; m++) {
k = binbig[ibin][m];
big = &biglist[k];
j = big->index;
+ if (j == jlast) continue;
type = big->type;
- if (type == SPHERE) inside = inside_sphere(x[i],x[j],big);
- else if (type == ELLIPSOID) inside = inside_ellipsoid(x[i],x[j],big);
- else inside = inside_wall(x[i],j);
+ if (type == SPHERE)
+ inside = inside_sphere(x[i],x[j],big);
+ else if (type == ELLIPSOID)
+ inside = inside_ellipsoid(x[i],x[j],big);
+ else if (type == LINE)
+ inside = inside_line(x[i],x[j],v[i],v[j],big,dt);
+ else if (type == TRIANGLE)
+ inside = inside_tri(x[i],x[j],v[i],v[j],big,dt);
+ else
+ inside = inside_wall(x[i],j);
if (inside) {
if (type == SPHERE)
@@ -1275,6 +1321,12 @@ void FixSRD::collisions_multi()
else if (type == ELLIPSOID)
t_remain = collision_ellipsoid_exact(x[i],x[j],v[i],v[j],big,
xscoll,xbcoll,norm);
+ else if (type == LINE)
+ t_remain = collision_line_exact(x[i],x[j],v[i],v[j],big,dt,
+ xscoll,xbcoll,norm);
+ else if (type == TRIANGLE)
+ t_remain = collision_tri_exact(x[i],x[j],v[i],v[j],big,dt,
+ xscoll,xbcoll,norm);
else
t_remain = collision_wall_exact(x[i],j,v[i],xscoll,xbcoll,norm);
@@ -1318,7 +1370,7 @@ void FixSRD::collisions_multi()
}
if (t_first == 0.0) break;
- j = jfirst;
+ j = jlast = jfirst;
type = typefirst;
xscoll[0] = xscollfirst[0];
xscoll[1] = xscollfirst[1];
@@ -1331,17 +1383,11 @@ void FixSRD::collisions_multi()
norm[2] = normfirst[2];
if (collidestyle == SLIP) {
- if (type == SPHERE)
- slip_sphere(v[i],v[j],norm,vsnew);
- else if (type == ELLIPSOID)
- slip_ellipsoid(v[i],v[j],x[j],big,xscoll,norm,vsnew);
- else
- slip_wall(v[i],j,norm,vsnew);
+ if (type != WALL) slip(v[i],v[j],x[j],big,xscoll,norm,vsnew);
+ else slip_wall(v[i],j,norm,vsnew);
} else {
- if (type != WALL)
- noslip(v[i],v[j],x[j],big,xscoll,norm,vsnew);
- else
- noslip_wall(v[i],j,xscoll,norm,vsnew);
+ if (type != WALL) noslip(v[i],v[j],x[j],big,-1,xscoll,norm,vsnew);
+ else noslip(v[i],NULL,x[j],big,j,xscoll,norm,vsnew);
}
if (dimension == 2) vsnew[2] = 0.0;
@@ -1420,6 +1466,276 @@ int FixSRD::inside_ellipsoid(double *xs, double *xb, Big *big)
return 0;
}
+/* ----------------------------------------------------------------------
+ check if SRD particle S is inside line big particle B
+ collision only possible if:
+ S starts on positive side of infinite line,
+ which means it will collide with outside of rigid body made of lines
+ since line segments have outward normals,
+ when vector from first to last point is crossed with +z axis
+ S ends on negative side of infinite line
+ unlike most other inside() routines, then calculate exact collision:
+ solve for collision pt along infinite line
+ collision if pt is within endpoints of B
+------------------------------------------------------------------------- */
+
+int FixSRD::inside_line(double *xs, double *xb, double *vs, double *vb,
+ Big *big, double dt_step)
+{
+ double pmc0[2],pmc1[2],n0[2],n1[2];
+ double n1_n0[2],pmc1_pmc0[2];
+
+ // 1 and 2 = start and end of timestep
+ // pmc = P - C, where P = position of S, C = position of B
+ // n = normal to line = [-sin(theta),cos(theta)], theta = orientation of B
+ // (P-C) dot N = side of line that S is on
+ // side0 = -1,0,1 for which side of line B that S is on at start of step
+ // side1 = -1,0,1 for which side of line B that S is on at end of step
+
+ xs1[0] = xs[0];
+ xs1[1] = xs[1];
+ xb1[0] = xb[0];
+ xb1[1] = xb[1];
+
+ xs0[0] = xs1[0] - dt_step*vs[0];
+ xs0[1] = xs1[1] - dt_step*vs[1];
+ xb0[0] = xb1[0] - dt_step*vb[0];
+ xb0[1] = xb1[1] - dt_step*vb[1];
+
+ theta1 = big->theta;
+ theta0 = theta1 - dt_step*big->omega[2];
+
+ pmc0[0] = xs0[0] - xb0[0];
+ pmc0[1] = xs0[1] - xb0[1];
+ n0[0] = sin(theta0);
+ n0[1] = -cos(theta0);
+
+ pmc1[0] = xs1[0] - xb1[0];
+ pmc1[1] = xs1[1] - xb1[1];
+ n1[0] = sin(theta1);
+ n1[1] = -cos(theta1);
+
+ double side0 = pmc0[0]*n0[0] + pmc0[1]*n0[1];
+ double side1 = pmc1[0]*n1[0] + pmc1[1]*n1[1];
+
+ if (side0 <= 0.0 || side1 >= 0.0) return 0;
+
+ // solve for time t (0 to 1) at which moving particle
+ // crosses infinite moving/rotating line
+
+ // Newton-Raphson solve of full non-linear parametric equation
+
+ tfraction = newton_raphson(0.0,1.0);
+
+ // quadratic equation solve of approximate parametric equation
+
+ /*
+ n1_n0[0] = n1[0]-n0[0]; n1_n0[1] = n1[1]-n0[1];
+ pmc1_pmc0[0] = pmc1[0]-pmc0[0]; pmc1_pmc0[1] = pmc1[1]-pmc0[1];
+
+ double a = pmc1_pmc0[0]*n1_n0[0] + pmc1_pmc0[1]*n1_n0[1];
+ double b = pmc1_pmc0[0]*n0[0] + pmc1_pmc0[1]*n0[1] +
+ n1_n0[0]*pmc0[0] + n1_n0[1]*pmc0[1];
+ double c = pmc0[0]*n0[0] + pmc0[1]*n0[1];
+
+ if (a == 0.0) {
+ double dot0 = pmc0[0]*n0[0] + pmc0[1]*n0[1];
+ double dot1 = pmc1[0]*n0[0] + pmc1[1]*n0[1];
+ double root = -dot0 / (dot1 - dot0);
+ //printf("Linear root: %g %g\n",root,tfraction);
+ tfraction = root;
+
+ } else {
+
+ double term = sqrt(b*b - 4.0*a*c);
+ double root1 = (-b + term) / (2.0*a);
+ double root2 = (-b - term) / (2.0*a);
+
+ //printf("ABC vecs: %g %g: %g %g\n",
+ // pmc1_pmc0[0],pmc1_pmc0[1],n1_n0[0],n1_n0[1]);
+ //printf("ABC vecs: %g %g: %g %g: %g %g %g\n",
+ // n0[0],n0[1],n1[0],n1[1],theta0,theta1,big->omega[2]);
+ //printf("ABC root: %g %g %g: %g %g %g\n",a,b,c,root1,root2,tfraction);
+
+ if (0.0 <= root1 && root1 <= 1.0) tfraction = root1;
+ else if (0.0 <= root2 && root2 <= 1.0) tfraction = root2;
+ else error->one(FLERR,"Bad quadratic solve for particle/line collision");
+ }
+ */
+
+ // check if collision pt is within line segment at collision time
+
+ xsc[0] = xs0[0] + tfraction*(xs1[0]-xs0[0]);
+ xsc[1] = xs0[1] + tfraction*(xs1[1]-xs0[1]);
+ xbc[0] = xb0[0] + tfraction*(xb1[0]-xb0[0]);
+ xbc[1] = xb0[1] + tfraction*(xb1[1]-xb0[1]);
+ double delx = xsc[0] - xbc[0];
+ double dely = xsc[1] - xbc[1];
+ double rsq = delx*delx + dely*dely;
+ if (rsq > 0.25*big->length*big->length) return 0;
+
+ //nbc[0] = n0[0] + tfraction*(n1[0]-n0[0]);
+ //nbc[1] = n0[1] + tfraction*(n1[1]-n0[1]);
+
+ nbc[0] = sin(theta0 + tfraction*(theta1-theta0));
+ nbc[1] = -cos(theta0 + tfraction*(theta1-theta0));
+
+ return 1;
+}
+
+/* ----------------------------------------------------------------------
+ check if SRD particle S is inside triangle big particle B
+ collision only possible if:
+ S starts on positive side of triangle plane,
+ which means it will collide with outside of rigid body made of tris
+ since triangles have outward normals,
+ S ends on negative side of triangle plane
+ unlike most other inside() routines, then calculate exact collision:
+ solve for collision pt on triangle plane
+ collision if pt is inside triangle B
+------------------------------------------------------------------------- */
+
+int FixSRD::inside_tri(double *xs, double *xb, double *vs, double *vb,
+ Big *big, double dt_step)
+{
+ double pmc0[3],pmc1[3],n0[3];
+ double n1_n0[3],pmc1_pmc0[3];
+ double excoll[3],eycoll[3],ezcoll[3];
+ double dc1[3],dc2[3],dc3[3];
+ double c1[3],c2[3],c3[3];
+ double c2mc1[3],c3mc2[3],c1mc3[3];
+ double pvec[3],xproduct[3];
+
+ // 1 and 2 = start and end of timestep
+ // pmc = P - C, where P = position of S, C = position of B
+ // n = normal to triangle
+ // (P-C) dot N = side of tri that S is on
+ // side0 = -1,0,1 for which side of tri B that S is on at start of step
+ // side1 = -1,0,1 for which side of tri B that S is on at end of step
+
+ double *omega = big->omega;
+ double *n1 = big->norm;
+
+ n0[0] = n1[0] - dt_step*(omega[1]*n1[2] - omega[2]*n1[1]);
+ n0[1] = n1[1] - dt_step*(omega[2]*n1[0] - omega[0]*n1[2]);
+ n0[2] = n1[2] - dt_step*(omega[0]*n1[1] - omega[1]*n1[0]);
+
+ pmc0[0] = xs[0] - dt_step*vs[0] - xb[0] + dt_step*vb[0];
+ pmc0[1] = xs[1] - dt_step*vs[1] - xb[1] + dt_step*vb[1];
+ pmc0[2] = xs[2] - dt_step*vs[2] - xb[2] + dt_step*vb[2];
+ pmc1[0] = xs[0] - xb[0];
+ pmc1[1] = xs[1] - xb[1];
+ pmc1[2] = xs[2] - xb[2];
+
+ double side0 = MathExtra::dot3(pmc0,n0);
+ double side1 = MathExtra::dot3(pmc1,n1);
+
+ if (side0 <= 0.0 || side1 >= 0.0) return 0;
+
+ // solve for time t (0 to 1) at which moving particle
+ // crosses moving/rotating tri
+ // quadratic equation solve of approximate parametric equation
+
+ n1_n0[0] = n1[0]-n0[0];
+ n1_n0[1] = n1[1]-n0[1];
+ n1_n0[2] = n1[2]-n0[2];
+ pmc1_pmc0[0] = pmc1[0]-pmc0[0];
+ pmc1_pmc0[1] = pmc1[1]-pmc0[1];
+ pmc1_pmc0[2] = pmc1[2]-pmc0[2];
+
+ double a = MathExtra::dot3(pmc1_pmc0,n1_n0);
+ double b = MathExtra::dot3(pmc1_pmc0,n0) + MathExtra::dot3(n1_n0,pmc0);
+ double c = MathExtra::dot3(pmc0,n0);
+
+ if (a == 0.0) {
+ double dot0 = MathExtra::dot3(pmc0,n0);
+ double dot1 = MathExtra::dot3(pmc1,n0);
+ double root = -dot0 / (dot1 - dot0);
+ tfraction = root;
+ } else {
+ double term = sqrt(b*b - 4.0*a*c);
+ double root1 = (-b + term) / (2.0*a);
+ double root2 = (-b - term) / (2.0*a);
+ if (0.0 <= root1 && root1 <= 1.0) tfraction = root1;
+ else if (0.0 <= root2 && root2 <= 1.0) tfraction = root2;
+ else error->one(FLERR,"Bad quadratic solve for particle/tri collision");
+ }
+
+ // calculate position/orientation of S and B at collision time
+ // dt = time previous to now at which collision occurs
+ // point = S position in plane of triangle at collision time
+ // Excoll,Eycoll,Ezcoll = orientation of tri at collision time
+ // c1,c2,c3 = corner points of triangle at collision time
+ // nbc = normal to plane of triangle at collision time
+
+ AtomVecTri::Bonus *tbonus;
+ tbonus = avec_tri->bonus;
+
+ double *ex = big->ex;
+ double *ey = big->ey;
+ double *ez = big->ez;
+ int index = atom->tri[big->index];
+ double *c1body = tbonus[index].c1;
+ double *c2body = tbonus[index].c2;
+ double *c3body = tbonus[index].c3;
+
+ double dt = (1.0-tfraction)*dt_step;
+
+ xsc[0] = xs[0] - dt*vs[0];
+ xsc[1] = xs[1] - dt*vs[1];
+ xsc[2] = xs[2] - dt*vs[2];
+ xbc[0] = xb[0] - dt*vb[0];
+ xbc[1] = xb[1] - dt*vb[1];
+ xbc[2] = xb[2] - dt*vb[2];
+
+ excoll[0] = ex[0] - dt*(omega[1]*ex[2] - omega[2]*ex[1]);
+ excoll[1] = ex[1] - dt*(omega[2]*ex[0] - omega[0]*ex[2]);
+ excoll[2] = ex[2] - dt*(omega[0]*ex[1] - omega[1]*ex[0]);
+
+ eycoll[0] = ey[0] - dt*(omega[1]*ey[2] - omega[2]*ey[1]);
+ eycoll[1] = ey[1] - dt*(omega[2]*ey[0] - omega[0]*ey[2]);
+ eycoll[2] = ey[2] - dt*(omega[0]*ey[1] - omega[1]*ey[0]);
+
+ ezcoll[0] = ez[0] - dt*(omega[1]*ez[2] - omega[2]*ez[1]);
+ ezcoll[1] = ez[1] - dt*(omega[2]*ez[0] - omega[0]*ez[2]);
+ ezcoll[2] = ez[2] - dt*(omega[0]*ez[1] - omega[1]*ez[0]);
+
+ MathExtra::matvec(excoll,eycoll,ezcoll,c1body,dc1);
+ MathExtra::matvec(excoll,eycoll,ezcoll,c2body,dc2);
+ MathExtra::matvec(excoll,eycoll,ezcoll,c3body,dc3);
+
+ MathExtra::add3(xbc,dc1,c1);
+ MathExtra::add3(xbc,dc2,c2);
+ MathExtra::add3(xbc,dc3,c3);
+
+ MathExtra::sub3(c2,c1,c2mc1);
+ MathExtra::sub3(c3,c2,c3mc2);
+ MathExtra::sub3(c1,c3,c1mc3);
+
+ MathExtra::cross3(c2mc1,c3mc2,nbc);
+ MathExtra::norm3(nbc);
+
+ // check if collision pt is within triangle
+ // pvec = vector from tri vertex to intersection point
+ // xproduct = cross product of edge vec with pvec
+ // if dot product of xproduct with nbc < 0.0 for any of 3 edges,
+ // intersection point is outside tri
+
+ MathExtra::sub3(xsc,c1,pvec);
+ MathExtra::cross3(c2mc1,pvec,xproduct);
+ if (MathExtra::dot3(xproduct,nbc) < 0.0) return 0;
+
+ MathExtra::sub3(xsc,c2,pvec);
+ MathExtra::cross3(c3mc2,pvec,xproduct);
+ if (MathExtra::dot3(xproduct,nbc) < 0.0) return 0;
+
+ MathExtra::sub3(xsc,c3,pvec);
+ MathExtra::cross3(c1mc3,pvec,xproduct);
+ if (MathExtra::dot3(xproduct,nbc) < 0.0) return 0;
+
+ return 1;
+}
+
/* ----------------------------------------------------------------------
check if SRD particle S is inside wall IWALL
------------------------------------------------------------------------- */
@@ -1451,7 +1767,7 @@ double FixSRD::collision_sphere_exact(double *xs, double *xb,
double vs_dot_vs,vb_dot_vb,vs_dot_vb;
double vs_dot_xb,vb_dot_xs,vs_dot_xs,vb_dot_xb;
double xs_dot_xs,xb_dot_xb,xs_dot_xb;
- double a,b,c,scale,dt;
+ double a,b,c,scale;
vs_dot_vs = vs[0]*vs[0] + vs[1]*vs[1] + vs[2]*vs[2];
vb_dot_vb = vb[0]*vb[0] + vb[1]*vb[1] + vb[2]*vb[2];
@@ -1470,7 +1786,7 @@ double FixSRD::collision_sphere_exact(double *xs, double *xb,
b = 2.0 * (vs_dot_xb + vb_dot_xs - vs_dot_xs - vb_dot_xb);
c = xs_dot_xs + xb_dot_xb - 2.0*xs_dot_xb - big->radsq;
- dt = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*a);
+ double dt = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*a);
xscoll[0] = xs[0] - dt*vs[0];
xscoll[1] = xs[1] - dt*vs[1];
@@ -1541,7 +1857,7 @@ double FixSRD::collision_ellipsoid_exact(double *xs, double *xb,
double vs_vb[3],xs_xb[3],omega_ex[3],omega_ey[3],omega_ez[3];
double excoll[3],eycoll[3],ezcoll[3],delta[3],xbody[3],nbody[3];
double ax,bx,cx,ay,by,cy,az,bz,cz;
- double a,b,c,dt,scale;
+ double a,b,c,scale;
double *omega = big->omega;
double *ex = big->ex;
@@ -1551,17 +1867,9 @@ double FixSRD::collision_ellipsoid_exact(double *xs, double *xb,
vs_vb[0] = vs[0]-vb[0]; vs_vb[1] = vs[1]-vb[1]; vs_vb[2] = vs[2]-vb[2];
xs_xb[0] = xs[0]-xb[0]; xs_xb[1] = xs[1]-xb[1]; xs_xb[2] = xs[2]-xb[2];
- omega_ex[0] = omega[1]*ex[2] - omega[2]*ex[1];
- omega_ex[1] = omega[2]*ex[0] - omega[0]*ex[2];
- omega_ex[2] = omega[0]*ex[1] - omega[1]*ex[0];
-
- omega_ey[0] = omega[1]*ey[2] - omega[2]*ey[1];
- omega_ey[1] = omega[2]*ey[0] - omega[0]*ey[2];
- omega_ey[2] = omega[0]*ey[1] - omega[1]*ey[0];
-
- omega_ez[0] = omega[1]*ez[2] - omega[2]*ez[1];
- omega_ez[1] = omega[2]*ez[0] - omega[0]*ez[2];
- omega_ez[2] = omega[0]*ez[1] - omega[1]*ez[0];
+ MathExtra::cross3(omega,ex,omega_ex);
+ MathExtra::cross3(omega,ey,omega_ey);
+ MathExtra::cross3(omega,ez,omega_ez);
ax = vs_vb[0]*omega_ex[0] + vs_vb[1]*omega_ex[1] + vs_vb[2]*omega_ex[2];
bx = -(vs_vb[0]*ex[0] + vs_vb[1]*ex[1] + vs_vb[2]*ex[2]);
@@ -1586,7 +1894,7 @@ double FixSRD::collision_ellipsoid_exact(double *xs, double *xb,
c = cx*cx*big->aradsqinv + cy*cy*big->bradsqinv +
cz*cz*big->cradsqinv - 1.0;
- dt = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*a);
+ double dt = (-b + sqrt(b*b - 4.0*a*c)) / (2.0*a);
xscoll[0] = xs[0] - dt*vs[0];
xscoll[1] = xs[1] - dt*vs[1];
@@ -1602,37 +1910,27 @@ double FixSRD::collision_ellipsoid_exact(double *xs, double *xb,
// norm = normal in space frame
// only worry about normalizing final norm vector
- excoll[0] = ex[0] - dt * (omega[1]*ex[2] - omega[2]*ex[1]);
- excoll[1] = ex[1] - dt * (omega[2]*ex[0] - omega[0]*ex[2]);
- excoll[2] = ex[2] - dt * (omega[0]*ex[1] - omega[1]*ex[0]);
-
- eycoll[0] = ey[0] - dt * (omega[1]*ey[2] - omega[2]*ey[1]);
- eycoll[1] = ey[1] - dt * (omega[2]*ey[0] - omega[0]*ey[2]);
- eycoll[2] = ey[2] - dt * (omega[0]*ey[1] - omega[1]*ey[0]);
+ excoll[0] = ex[0] - dt*(omega[1]*ex[2] - omega[2]*ex[1]);
+ excoll[1] = ex[1] - dt*(omega[2]*ex[0] - omega[0]*ex[2]);
+ excoll[2] = ex[2] - dt*(omega[0]*ex[1] - omega[1]*ex[0]);
- ezcoll[0] = ez[0] - dt * (omega[1]*ez[2] - omega[2]*ez[1]);
- ezcoll[1] = ez[1] - dt * (omega[2]*ez[0] - omega[0]*ez[2]);
- ezcoll[2] = ez[2] - dt * (omega[0]*ez[1] - omega[1]*ez[0]);
+ eycoll[0] = ey[0] - dt*(omega[1]*ey[2] - omega[2]*ey[1]);
+ eycoll[1] = ey[1] - dt*(omega[2]*ey[0] - omega[0]*ey[2]);
+ eycoll[2] = ey[2] - dt*(omega[0]*ey[1] - omega[1]*ey[0]);
- delta[0] = xscoll[0] - xbcoll[0];
- delta[1] = xscoll[1] - xbcoll[1];
- delta[2] = xscoll[2] - xbcoll[2];
+ ezcoll[0] = ez[0] - dt*(omega[1]*ez[2] - omega[2]*ez[1]);
+ ezcoll[1] = ez[1] - dt*(omega[2]*ez[0] - omega[0]*ez[2]);
+ ezcoll[2] = ez[2] - dt*(omega[0]*ez[1] - omega[1]*ez[0]);
- xbody[0] = delta[0]*excoll[0] + delta[1]*excoll[1] + delta[2]*excoll[2];
- xbody[1] = delta[0]*eycoll[0] + delta[1]*eycoll[1] + delta[2]*eycoll[2];
- xbody[2] = delta[0]*ezcoll[0] + delta[1]*ezcoll[1] + delta[2]*ezcoll[2];
+ MathExtra::sub3(xscoll,xbcoll,delta);
+ MathExtra::transpose_matvec(excoll,eycoll,ezcoll,delta,xbody);
nbody[0] = xbody[0]*big->aradsqinv;
nbody[1] = xbody[1]*big->bradsqinv;
nbody[2] = xbody[2]*big->cradsqinv;
- norm[0] = excoll[0]*nbody[0] + eycoll[0]*nbody[1] + ezcoll[0]*nbody[2];
- norm[1] = excoll[1]*nbody[0] + eycoll[1]*nbody[1] + ezcoll[1]*nbody[2];
- norm[2] = excoll[2]*nbody[0] + eycoll[2]*nbody[1] + ezcoll[2]*nbody[2];
- scale = 1.0/sqrt(norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]);
- norm[0] *= scale;
- norm[1] *= scale;
- norm[2] *= scale;
+ MathExtra::matvec(excoll,eycoll,ezcoll,nbody,norm);
+ MathExtra::norm3(norm);
return dt;
}
@@ -1640,6 +1938,10 @@ double FixSRD::collision_ellipsoid_exact(double *xs, double *xb,
/* ----------------------------------------------------------------------
collision of SRD particle S with surface of ellipsoidal big particle B
inexact because just push SRD to surface of big particle at end of step
+ time of collision = end of step
+ xscoll = collision pt = position of SRD at time of collision
+ xbcoll = xb = position of big particle at time of collision
+ norm = surface normal of collision pt at time of collision
------------------------------------------------------------------------- */
void FixSRD::collision_ellipsoid_inexact(double *xs, double *xb,
@@ -1648,22 +1950,18 @@ void FixSRD::collision_ellipsoid_inexact(double *xs, double *xb,
double *norm)
{
double xs_xb[3],delta[3],xbody[3],nbody[3];
- double x,y,z,scale;
double *ex = big->ex;
double *ey = big->ey;
double *ez = big->ez;
- xs_xb[0] = xs[0] - xb[0];
- xs_xb[1] = xs[1] - xb[1];
- xs_xb[2] = xs[2] - xb[2];
+ MathExtra::sub3(xs,xb,xs_xb);
+ double x = MathExtra::dot3(xs_xb,ex);
+ double y = MathExtra::dot3(xs_xb,ey);
+ double z = MathExtra::dot3(xs_xb,ez);
- x = xs_xb[0]*ex[0] + xs_xb[1]*ex[1] + xs_xb[2]*ex[2];
- y = xs_xb[0]*ey[0] + xs_xb[1]*ey[1] + xs_xb[2]*ey[2];
- z = xs_xb[0]*ez[0] + xs_xb[1]*ez[1] + xs_xb[2]*ez[2];
-
- scale = 1.0/sqrt(x*x*big->aradsqinv + y*y*big->bradsqinv +
- z*z*big->cradsqinv);
+ double scale = 1.0/sqrt(x*x*big->aradsqinv + y*y*big->bradsqinv +
+ z*z*big->cradsqinv);
x *= scale;
y *= scale;
z *= scale;
@@ -1681,25 +1979,73 @@ void FixSRD::collision_ellipsoid_inexact(double *xs, double *xb,
// norm = normal in space frame
// only worry about normalizing final norm vector
- delta[0] = xscoll[0] - xbcoll[0];
- delta[1] = xscoll[1] - xbcoll[1];
- delta[2] = xscoll[2] - xbcoll[2];
-
- xbody[0] = delta[0]*ex[0] + delta[1]*ex[1] + delta[2]*ex[2];
- xbody[1] = delta[0]*ey[0] + delta[1]*ey[1] + delta[2]*ey[2];
- xbody[2] = delta[0]*ez[0] + delta[1]*ez[1] + delta[2]*ez[2];
+ MathExtra::sub3(xscoll,xbcoll,delta);
+ MathExtra::transpose_matvec(ex,ey,ez,delta,xbody);
nbody[0] = xbody[0]*big->aradsqinv;
nbody[1] = xbody[1]*big->bradsqinv;
nbody[2] = xbody[2]*big->cradsqinv;
- norm[0] = ex[0]*nbody[0] + ey[0]*nbody[1] + ez[0]*nbody[2];
- norm[1] = ex[1]*nbody[0] + ey[1]*nbody[1] + ez[1]*nbody[2];
- norm[2] = ex[2]*nbody[0] + ey[2]*nbody[1] + ez[2]*nbody[2];
- scale = 1.0/sqrt(norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]);
- norm[0] *= scale;
- norm[1] *= scale;
- norm[2] *= scale;
+ MathExtra::matvec(ex,ey,ez,nbody,norm);
+ MathExtra::norm3(norm);
+}
+
+/* ----------------------------------------------------------------------
+ collision of SRD particle S with surface of line big particle B
+ exact because compute time of collision
+ dt = time previous to now at which collision occurs
+ xscoll = collision pt = position of SRD at time of collision
+ xbcoll = position of big particle at time of collision
+ norm = surface normal of collision pt at time of collision
+------------------------------------------------------------------------- */
+
+double FixSRD::collision_line_exact(double *xs, double *xb,
+ double *vs, double *vb, Big *big,
+ double dt_step,
+ double *xscoll, double *xbcoll,
+ double *norm)
+{
+ xscoll[0] = xsc[0];
+ xscoll[1] = xsc[1];
+ xscoll[2] = 0.0;
+ xbcoll[0] = xbc[0];
+ xbcoll[1] = xbc[1];
+ xbcoll[2] = 0.0;
+
+ norm[0] = nbc[0];
+ norm[1] = nbc[1];
+ norm[2] = 0.0;
+
+ return (1.0-tfraction)*dt_step;
+}
+
+/* ----------------------------------------------------------------------
+ collision of SRD particle S with surface of triangle big particle B
+ exact because compute time of collision
+ dt = time previous to now at which collision occurs
+ xscoll = collision pt = position of SRD at time of collision
+ xbcoll = position of big particle at time of collision
+ norm = surface normal of collision pt at time of collision
+------------------------------------------------------------------------- */
+
+double FixSRD::collision_tri_exact(double *xs, double *xb,
+ double *vs, double *vb, Big *big,
+ double dt_step,
+ double *xscoll, double *xbcoll,
+ double *norm)
+{
+ xscoll[0] = xsc[0];
+ xscoll[1] = xsc[1];
+ xscoll[2] = xsc[2];
+ xbcoll[0] = xbc[0];
+ xbcoll[1] = xbc[1];
+ xbcoll[2] = xbc[2];
+
+ norm[0] = nbc[0];
+ norm[1] = nbc[1];
+ norm[2] = nbc[2];
+
+ return (1.0-tfraction)*dt_step;
}
/* ----------------------------------------------------------------------
@@ -1707,6 +2053,7 @@ void FixSRD::collision_ellipsoid_inexact(double *xs, double *xb,
exact because compute time of collision
dt = time previous to now at which collision occurs
xscoll = collision pt = position of SRD at time of collision
+ xbcoll = position of wall at time of collision
norm = surface normal of collision pt at time of collision
------------------------------------------------------------------------- */
@@ -1737,6 +2084,7 @@ double FixSRD::collision_wall_exact(double *xs, int iwall, double *vs,
inexact because just push SRD to surface of wall at end of step
time of collision = end of step
xscoll = collision pt = position of SRD at time of collision
+ xbcoll = position of wall at time of collision
norm = surface normal of collision pt at time of collision
------------------------------------------------------------------------- */
@@ -1760,54 +2108,18 @@ void FixSRD::collision_wall_inexact(double *xs, int iwall, double *xscoll,
}
/* ----------------------------------------------------------------------
- SLIP collision with sphere
- vs = velocity of SRD, vb = velocity of BIG
- norm = unit normal from surface of BIG at collision pt
- v of BIG particle in direction of surf normal is added to v of SRD
- return vsnew of SRD
-------------------------------------------------------------------------- */
-
-void FixSRD::slip_sphere(double *vs, double *vb, double *norm, double *vsnew)
-{
- double r1,r2,vnmag,vs_dot_n,vsurf_dot_n;
- double tangent[3];
-
- while (1) {
- r1 = sigma * random->gaussian();
- r2 = sigma * random->gaussian();
- vnmag = sqrt(r1*r1 + r2*r2);
- if (vnmag*vnmag <= vmaxsq) break;
- }
-
- vs_dot_n = vs[0]*norm[0] + vs[1]*norm[1] + vs[2]*norm[2];
-
- tangent[0] = vs[0] - vs_dot_n*norm[0];
- tangent[1] = vs[1] - vs_dot_n*norm[1];
- tangent[2] = vs[2] - vs_dot_n*norm[2];
-
- // vsurf = velocity of collision pt = translation/rotation of BIG particle
- // for sphere, only vb (translation) can contribute in normal direction
-
- vsurf_dot_n = vb[0]*norm[0] + vb[1]*norm[1] + vb[2]*norm[2];
-
- vsnew[0] = (vnmag+vsurf_dot_n)*norm[0] + tangent[0];
- vsnew[1] = (vnmag+vsurf_dot_n)*norm[1] + tangent[1];
- vsnew[2] = (vnmag+vsurf_dot_n)*norm[2] + tangent[2];
-}
-
-/* ----------------------------------------------------------------------
- SLIP collision with ellipsoid
+ SLIP collision with BIG particle with omega
vs = velocity of SRD, vb = velocity of BIG
xb = position of BIG, omega = rotation of BIG
xsurf = collision pt on surf of BIG
norm = unit normal from surface of BIG at collision pt
v of BIG particle in direction of surf normal is added to v of SRD
- includes component due to rotation of ellipsoid
+ includes component due to rotation of BIG
return vsnew of SRD
------------------------------------------------------------------------- */
-void FixSRD::slip_ellipsoid(double *vs, double *vb, double *xb, Big *big,
- double *xsurf, double *norm, double *vsnew)
+void FixSRD::slip(double *vs, double *vb, double *xb, Big *big,
+ double *xsurf, double *norm, double *vsnew)
{
double r1,r2,vnmag,vs_dot_n,vsurf_dot_n;
double tangent[3],vsurf[3];
@@ -1827,6 +2139,8 @@ void FixSRD::slip_ellipsoid(double *vs, double *vb, double *xb, Big *big,
tangent[2] = vs[2] - vs_dot_n*norm[2];
// vsurf = velocity of collision pt = translation/rotation of BIG particle
+ // NOTE: for sphere could just use vsurf = vb, since w x (xsurf-xb)
+ // is orthogonal to norm and thus doesn't contribute to vsurf_dot_n
vsurf[0] = vb[0] + omega[1]*(xsurf[2]-xb[2]) - omega[2]*(xsurf[1]-xb[1]);
vsurf[1] = vb[1] + omega[2]*(xsurf[0]-xb[0]) - omega[0]*(xsurf[2]-xb[2]);
@@ -1887,7 +2201,7 @@ void FixSRD::slip_wall(double *vs, int iwall, double *norm, double *vsnew)
}
/* ----------------------------------------------------------------------
- NO-SLIP collision with sphere or ellipsoid
+ NO-SLIP collision with BIG particle including WALL
vs = velocity of SRD, vb = velocity of BIG
xb = position of BIG, omega = rotation of BIG
xsurf = collision pt on surf of BIG
@@ -1897,12 +2211,11 @@ void FixSRD::slip_wall(double *vs, int iwall, double *norm, double *vsnew)
return vsnew of SRD
------------------------------------------------------------------------- */
-void FixSRD::noslip(double *vs, double *vb, double *xb, Big *big,
+void FixSRD::noslip(double *vs, double *vb, double *xb, Big *big, int iwall,
double *xsurf, double *norm, double *vsnew)
{
double vs_dot_n,scale,r1,r2,vnmag,vtmag1,vtmag2;
double tangent1[3],tangent2[3];
- double *omega = big->omega;
vs_dot_n = vs[0]*norm[0] + vs[1]*norm[1] + vs[2]*norm[2];
@@ -1932,60 +2245,20 @@ void FixSRD::noslip(double *vs, double *vb, double *xb, Big *big,
vsnew[1] = vnmag*norm[1] + vtmag1*tangent1[1] + vtmag2*tangent2[1];
vsnew[2] = vnmag*norm[2] + vtmag1*tangent1[2] + vtmag2*tangent2[2];
- // add in velocity of collision pt = translation/rotation of BIG particle
+ // add in velocity of collision pt
+ // for WALL: velocity of wall in one dim
+ // else: translation/rotation of BIG particle
- vsnew[0] += vb[0] + omega[1]*(xsurf[2]-xb[2]) - omega[2]*(xsurf[1]-xb[1]);
- vsnew[1] += vb[1] + omega[2]*(xsurf[0]-xb[0]) - omega[0]*(xsurf[2]-xb[2]);
- vsnew[2] += vb[2] + omega[0]*(xsurf[1]-xb[1]) - omega[1]*(xsurf[0]-xb[0]);
-}
+ if (big->type == WALL) {
+ int dim = wallwhich[iwall] / 2;
+ vsnew[dim] += vwall[iwall];
-/* ----------------------------------------------------------------------
- NO-SLIP collision with wall IWALL
- vs = velocity of SRD
- xsurf = collision pt on surf of WALL
- norm = unit normal from WALL at collision pt
- v of collision pt is added to v of SRD
- return vsnew of SRD
-------------------------------------------------------------------------- */
-
-void FixSRD::noslip_wall(double *vs, int iwall,
- double *xsurf, double *norm, double *vsnew)
-{
- double vs_dot_n,scale,r1,r2,vnmag,vtmag1,vtmag2;
- double tangent1[3],tangent2[3];
-
- vs_dot_n = vs[0]*norm[0] + vs[1]*norm[1] + vs[2]*norm[2];
-
- tangent1[0] = vs[0] - vs_dot_n*norm[0];
- tangent1[1] = vs[1] - vs_dot_n*norm[1];
- tangent1[2] = vs[2] - vs_dot_n*norm[2];
- scale = 1.0/sqrt(tangent1[0]*tangent1[0] + tangent1[1]*tangent1[1] +
- tangent1[2]*tangent1[2]);
- tangent1[0] *= scale;
- tangent1[1] *= scale;
- tangent1[2] *= scale;
-
- tangent2[0] = norm[1]*tangent1[2] - norm[2]*tangent1[1];
- tangent2[1] = norm[2]*tangent1[0] - norm[0]*tangent1[2];
- tangent2[2] = norm[0]*tangent1[1] - norm[1]*tangent1[0];
-
- while (1) {
- r1 = sigma * random->gaussian();
- r2 = sigma * random->gaussian();
- vnmag = sqrt(r1*r1 + r2*r2);
- vtmag1 = sigma * random->gaussian();
- vtmag2 = sigma * random->gaussian();
- if (vnmag*vnmag + vtmag1*vtmag1 + vtmag2*vtmag2 <= vmaxsq) break;
+ } else {
+ double *omega = big->omega;
+ vsnew[0] += vb[0] + omega[1]*(xsurf[2]-xb[2]) - omega[2]*(xsurf[1]-xb[1]);
+ vsnew[1] += vb[1] + omega[2]*(xsurf[0]-xb[0]) - omega[0]*(xsurf[2]-xb[2]);
+ vsnew[2] += vb[2] + omega[0]*(xsurf[1]-xb[1]) - omega[1]*(xsurf[0]-xb[0]);
}
-
- vsnew[0] = vnmag*norm[0] + vtmag1*tangent1[0] + vtmag2*tangent2[0];
- vsnew[1] = vnmag*norm[1] + vtmag1*tangent1[1] + vtmag2*tangent2[1];
- vsnew[2] = vnmag*norm[2] + vtmag1*tangent1[2] + vtmag2*tangent2[2];
-
- // add in velocity of collision pt = velocity of wall
-
- int dim = wallwhich[iwall] / 2;
- vsnew[dim] += vwall[iwall];
}
/* ----------------------------------------------------------------------
@@ -2027,6 +2300,7 @@ void FixSRD::force_torque(double *vsold, double *vsnew,
------------------------------------------------------------------------- */
void FixSRD::force_wall(double *vsold, double *vsnew, int iwall)
+
{
double dpdt[3];
@@ -2100,12 +2374,20 @@ void FixSRD::parameterize()
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
double *radius = atom->radius;
int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
int *mask = atom->mask;
int nlocal = atom->nlocal;
- any_ellipsoids = 0;
+ int any_ellipsoids = 0;
+ int any_lines = 0;
+ int any_tris = 0;
maxbigdiam = 0.0;
minbigdiam = BIG;
@@ -2123,6 +2405,20 @@ void FixSRD::parameterize()
minbigdiam = MIN(minbigdiam,2.0*shape[0]);
minbigdiam = MIN(minbigdiam,2.0*shape[1]);
minbigdiam = MIN(minbigdiam,2.0*shape[2]);
+ } else if (line && line[i] >= 0) {
+ any_lines = 1;
+ double length = lbonus[line[i]].length;
+ maxbigdiam = MAX(maxbigdiam,length);
+ minbigdiam = MIN(minbigdiam,length);
+ } else if (tri && tri[i] >= 0) {
+ any_tris = 1;
+ double length1 = MathExtra::len3(tbonus[tri[i]].c1);
+ double length2 = MathExtra::len3(tbonus[tri[i]].c2);
+ double length3 = MathExtra::len3(tbonus[tri[i]].c3);
+ double length = MAX(length1,length2);
+ length = MAX(length,length3);
+ maxbigdiam = MAX(maxbigdiam,length);
+ minbigdiam = MIN(minbigdiam,length);
} else
error->one(FLERR,"Big particle in fix srd cannot be point particle");
}
@@ -2137,10 +2433,20 @@ void FixSRD::parameterize()
int itmp = any_ellipsoids;
MPI_Allreduce(&itmp,&any_ellipsoids,1,MPI_INT,MPI_MAX,world);
+ itmp = any_lines;
+ MPI_Allreduce(&itmp,&any_lines,1,MPI_INT,MPI_MAX,world);
+ itmp = any_tris;
+ MPI_Allreduce(&itmp,&any_tris,1,MPI_INT,MPI_MAX,world);
+
+ if (any_lines && overlap == 0)
+ error->all(FLERR,"Cannot use lines with fix srd unless overlap is set");
+ if (any_tris && overlap == 0)
+ error->all(FLERR,"Cannot use tris with fix srd unless overlap is set");
- // big particles are only torqued if are ellipsoids or NOSLIP collisions
+ // big particles are only torqued if ellipsoids/lines/tris or NOSLIP
- if (any_ellipsoids == 0 && collidestyle == SLIP) torqueflag = 0;
+ if (any_ellipsoids == 0 && any_lines == 0 && any_tris == 0 &&
+ collidestyle == SLIP) torqueflag = 0;
else torqueflag = 1;
// mass of SRD particles, require monodispersity
@@ -2186,28 +2492,46 @@ void FixSRD::parameterize()
vmaxsq = vmax*vmax;
// volbig = total volume of all big particles
+ // LINE/TRI particles have no volume
+ // incorrect if part of rigid particles, so add fudge factor with WIDTH
// apply radfactor to reduce final volume
double volbig = 0.0;
+ double WIDTH = 1.0;
if (dimension == 3) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & biggroupbit) {
- if (radius && radius[i] > 0.0)
- volbig += 4.0/3.0*MY_PI*radius[i]*radius[i]*radius[i];
- else if (ellipsoid && ellipsoid[i] >= 0) {
+ if (radius && radius[i] > 0.0) {
+ double r = radfactor * radius[i];
+ volbig += 4.0/3.0*MY_PI * r*r*r;;
+ } else if (ellipsoid && ellipsoid[i] >= 0) {
double *shape = ebonus[ellipsoid[i]].shape;
- volbig += 4.0/3.0*MY_PI * shape[0]*shape[1]*shape[2];
+ volbig += 4.0/3.0*MY_PI * shape[0]*shape[1]*shape[2] *
+ radfactor*radfactor*radfactor;
+ } else if (tri && tri[i] >= 0) {
+ double *c1 = tbonus[tri[i]].c1;
+ double *c2 = tbonus[tri[i]].c2;
+ double *c3 = tbonus[tri[i]].c3;
+ double c2mc1[3],c3mc1[3],cross[3];
+ MathExtra::sub3(c2,c1,c2mc1);
+ MathExtra::sub3(c3,c1,c3mc1);
+ MathExtra::cross3(c2mc1,c3mc1,cross);
+ volbig += 0.5 * MathExtra::len3(cross);
}
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & biggroupbit) {
- if (radius && radius[i] > 0.0)
- volbig += MY_PI*radius[i]*radius[i];
- else if (ellipsoid && ellipsoid[i] >= 0) {
+ if (radius && radius[i] > 0.0) {
+ double r = radfactor * radius[i];
+ volbig += MY_PI * r*r;
+ } else if (ellipsoid && ellipsoid[i] >= 0) {
double *shape = ebonus[ellipsoid[i]].shape;
- volbig += MY_PI*shape[0]*shape[1];
+ volbig += MY_PI * shape[0]*shape[1] * radfactor*radfactor;
+ } else if (line && line[i] >= 0) {
+ double length = lbonus[line[i]].length;
+ volbig += length * WIDTH;
}
}
}
@@ -2215,9 +2539,6 @@ void FixSRD::parameterize()
tmp = volbig;
MPI_Allreduce(&tmp,&volbig,1,MPI_DOUBLE,MPI_SUM,world);
- if (dimension == 3) volbig *= radfactor*radfactor*radfactor;
- else volbig *= radfactor*radfactor;
-
// particle counts
bigint mbig = 0;
@@ -2228,7 +2549,10 @@ void FixSRD::parameterize()
mass_big = 0.0;
for (int i = 0; i < nlocal; i++)
- if (mask[i] & biggroupbit) mass_big += rmass[i];
+ if (mask[i] & biggroupbit) {
+ if (rmass) mass_big += rmass[i];
+ else mass_big += mass[type[i]];
+ }
tmp = mass_big;
MPI_Allreduce(&tmp,&mass_big,1,MPI_DOUBLE,MPI_SUM,world);
@@ -2369,7 +2693,8 @@ void FixSRD::parameterize()
if (cubicflag == CUBIC_ERROR)
error->all(FLERR,"SRD bin size for fix srd differs from user request");
if (me == 0)
- error->warning(FLERR,"SRD bin size for fix srd differs from user request");
+ error->warning(FLERR,
+ "SRD bin size for fix srd differs from user request");
}
// error if lamda < 0.6 of SRD grid size and no shifting allowed
@@ -2400,30 +2725,46 @@ void FixSRD::parameterize()
/* ----------------------------------------------------------------------
set static parameters of each big particle, owned and ghost
called each reneighboring
- use radfactor in distance parameters
+ use radfactor in distance parameters as appropriate
------------------------------------------------------------------------- */
void FixSRD::big_static()
{
int i;
- double rad,arad,brad,crad;
- double *shape;
+ double rad,arad,brad,crad,length,length1,length2,length3;
+ double *shape,*c1,*c2,*c3;
+ double c2mc1[3],c3mc1[3];
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
double *radius = atom->radius;
int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
+ int *type = atom->type;
double skinhalf = 0.5 * neighbor->skin;
for (int k = 0; k < nbig; k++) {
i = biglist[k].index;
+
+ // sphere
+ // set radius and radsq and cutoff based on radius
+
if (radius && radius[i] > 0.0) {
biglist[k].type = SPHERE;
rad = radfactor*radius[i];
biglist[k].radius = rad;
biglist[k].radsq = rad*rad;
biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
+
+ // ellipsoid
+ // set abc radsqinv and cutoff based on max radius
+
} else if (ellipsoid && ellipsoid[i] >= 0) {
shape = ebonus[ellipsoid[i]].shape;
biglist[k].type = ELLIPSOID;
@@ -2436,33 +2777,67 @@ void FixSRD::big_static()
rad = MAX(arad,brad);
rad = MAX(rad,crad);
biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
+
+ // line
+ // set length and cutoff based on 1/2 length
+
+ } else if (line && line[i] >= 0) {
+ length = lbonus[line[i]].length;
+ biglist[k].type = LINE;
+ biglist[k].length = length;
+ rad = 0.5*length;
+ biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
+
+ // tri
+ // set normbody based on c1,c2,c3
+ // set cutoff based on point furthest from centroid
+
+ } else if (tri && tri[i] >= 0) {
+ biglist[k].type = TRIANGLE;
+ c1 = tbonus[tri[i]].c1;
+ c2 = tbonus[tri[i]].c2;
+ c3 = tbonus[tri[i]].c3;
+ MathExtra::sub3(c2,c1,c2mc1);
+ MathExtra::sub3(c3,c1,c3mc1);
+ MathExtra::cross3(c2mc1,c3mc1,biglist[k].normbody);
+ length1 = MathExtra::len3(c1);
+ length2 = MathExtra::len3(c2);
+ length3 = MathExtra::len3(c3);
+ rad = MAX(length1,length2);
+ rad = MAX(rad,length3);
+ biglist[k].cutbinsq = (rad+skinhalf) * (rad+skinhalf);
}
}
}
/* ----------------------------------------------------------------------
- set dynamics parameters of each big particle, owned and ghost
- for ELLIPSOID, need current omega and ex,ey,ez
- for SPHERE, need current omega from atom->omega or atom->angmom
+ set dynamic parameters of each big particle, owned and ghost
called each timestep
------------------------------------------------------------------------- */
void FixSRD::big_dynamic()
{
int i;
- double *shape,*quat;
+ double *shape,*quat,*inertia;
+ double inertiaone[3];
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
double **omega = atom->omega;
double **angmom = atom->angmom;
double *rmass = atom->rmass;
int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
for (int k = 0; k < nbig; k++) {
i = biglist[k].index;
- // sphere with omega
+ // sphere
// set omega from atom->omega directly
if (biglist[k].type == SPHERE) {
@@ -2470,15 +2845,45 @@ void FixSRD::big_dynamic()
biglist[k].omega[1] = omega[i][1];
biglist[k].omega[2] = omega[i][2];
- // ellipsoid with angmom
- // calculate ex,ey,ez and omega from quaternion and angmom
+ // ellipsoid
+ // set ex,ey,ez from quaternion
+ // set omega from angmom & ex,ey,ez
} else if (biglist[k].type == ELLIPSOID) {
- shape = ebonus[ellipsoid[i]].shape;
quat = ebonus[ellipsoid[i]].quat;
- exyz_from_q(quat,biglist[k].ex,biglist[k].ey,biglist[k].ez);
- omega_from_mq(angmom[i],biglist[k].ex,biglist[k].ey,biglist[k].ez,
- rmass[i],shape,biglist[k].omega);
+ MathExtra::q_to_exyz(quat,biglist[k].ex,biglist[k].ey,biglist[k].ez);
+ shape = ebonus[ellipsoid[i]].shape;
+ inertiaone[0] = EINERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);
+ inertiaone[1] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);
+ inertiaone[2] = EINERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);
+ MathExtra::angmom_to_omega(angmom[i],
+ biglist[k].ex,biglist[k].ey,biglist[k].ez,
+ inertiaone,biglist[k].omega);
+
+ // line
+ // set omega from atom->omega directly
+
+ } else if (biglist[k].type == LINE) {
+ biglist[k].theta = lbonus[line[i]].theta;
+ biglist[k].omega[0] = omega[i][0];
+ biglist[k].omega[1] = omega[i][1];
+ biglist[k].omega[2] = omega[i][2];
+
+ // tri
+ // set ex,ey,ez from quaternion
+ // set omega from angmom & ex,ey,ez
+ // set unit space-frame norm from body-frame norm
+
+ } else if (biglist[k].type == TRIANGLE) {
+ quat = tbonus[tri[i]].quat;
+ MathExtra::q_to_exyz(quat,biglist[k].ex,biglist[k].ey,biglist[k].ez);
+ inertia = tbonus[tri[i]].inertia;
+ MathExtra::angmom_to_omega(angmom[i],
+ biglist[k].ex,biglist[k].ey,biglist[k].ez,
+ inertia,biglist[k].omega);
+ MathExtra::matvec(biglist[k].ex,biglist[k].ey,biglist[k].ez,
+ biglist[k].normbody,biglist[k].norm);
+ MathExtra::norm3(biglist[k].norm);
}
}
}
@@ -2967,10 +3372,10 @@ void FixSRD::setup_velocity_shift(int ishift, int dynamic)
for (m = 0; m < nsend; m++) vbin[sendlist[m]].owner = 0;
}
- // if triclinic and streamflag:
- // set xctr (in lamda units) for all nbins so can compute bin vstream
+ // if tstat and deformflag:
+ // set xctr for all nbins in lamda units so can later compute vstream of bin
- if (triclinic && streamflag) {
+ if (tstat && deformflag) {
m = 0;
for (k = 0; k < nbinz; k++)
for (j = 0; j < nbiny; j++)
@@ -3182,55 +3587,6 @@ double FixSRD::bin_bin_distance(int i, int j, int k)
return (delx*delx + dely*dely + delz*delz);
}
-/* ----------------------------------------------------------------------
- compute space-frame ex,ey,ez from current quaternion q
- ex,ey,ez = space-frame coords of 1st,2nd,3rd principal axis
- operation is ex = q' d q = Q d, where d is (1,0,0) = 1st axis in body frame
-------------------------------------------------------------------------- */
-
-void FixSRD::exyz_from_q(double *q, double *ex, double *ey, double *ez)
-{
- ex[0] = q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3];
- ex[1] = 2.0 * (q[1]*q[2] + q[0]*q[3]);
- ex[2] = 2.0 * (q[1]*q[3] - q[0]*q[2]);
-
- ey[0] = 2.0 * (q[1]*q[2] - q[0]*q[3]);
- ey[1] = q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3];
- ey[2] = 2.0 * (q[2]*q[3] + q[0]*q[1]);
-
- ez[0] = 2.0 * (q[1]*q[3] + q[0]*q[2]);
- ez[1] = 2.0 * (q[2]*q[3] - q[0]*q[1]);
- ez[2] = q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3];
-}
-
-/* ----------------------------------------------------------------------
- compute omega from angular momentum
- w = omega = angular velocity in space frame
- wbody = angular velocity in body frame
- set wbody component to 0.0 if inertia component is 0.0
- otherwise body can spin easily around that axis
- project space-frame angular momentum onto body axes
- and divide by principal moments
-------------------------------------------------------------------------- */
-
-void FixSRD::omega_from_mq(double *m, double *ex, double *ey, double *ez,
- double mass, double *shape, double *w)
-{
- double inertia[3],wbody[3];
-
- inertia[0] = 0.2*mass * (shape[1]*shape[1]+shape[2]*shape[2]);
- inertia[1] = 0.2*mass * (shape[0]*shape[0]+shape[2]*shape[2]);
- inertia[2] = 0.2*mass * (shape[0]*shape[0]+shape[1]*shape[1]);
-
- wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / inertia[0];
- wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / inertia[1];
- wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / inertia[2];
-
- w[0] = wbody[0]*ex[0] + wbody[1]*ey[0] + wbody[2]*ez[0];
- w[1] = wbody[0]*ex[1] + wbody[1]*ey[1] + wbody[2]*ez[1];
- w[2] = wbody[0]*ex[2] + wbody[1]*ey[2] + wbody[2]*ez[2];
-}
-
/* ---------------------------------------------------------------------- */
int FixSRD::pack_reverse_comm(int n, int first, double *buf)
@@ -3365,6 +3721,90 @@ void FixSRD::velocity_stats(int groupnum)
/* ---------------------------------------------------------------------- */
+double FixSRD::newton_raphson(double t1, double t2)
+{
+ double f1,f2,df,tlo,thi;
+ lineside(t1,f1,df);
+ if (f1 < 0.0) {
+ tlo = t1;
+ thi = t2;
+ } else {
+ thi = t1;
+ tlo = t2;
+ }
+
+ double f;
+ double t = 0.5*(t1+t2);
+ double dtold = fabs(t2-t1);
+ double dt = dtold;
+ lineside(t,f,df);
+
+ double temp;
+ for (int i = 0; i < MAXITER; i++) {
+ if ((((t-thi)*df - f)*((t-tlo)*df - f) > 0.0) ||
+ (fabs(2.0*f) > fabs(dtold*df))) {
+ dtold = dt;
+ dt = 0.5 * (thi-tlo);
+ t = tlo + dt;
+ if (tlo == t) return t;
+ } else {
+ dtold = dt;
+ dt = f / df;
+ temp = t;
+ t -= dt;
+ if (temp == t) return t;
+ }
+ if (fabs(dt) < TOLERANCE) return t;
+ lineside(t,f,df);
+ if (f < 0.0) tlo = t;
+ else thi = t;
+ }
+
+ return t;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSRD::lineside(double t, double &f, double &df)
+{
+ double p[2],c[2];
+
+ p[0] = xs0[0] + (xs1[0]-xs0[0])*t;
+ p[1] = xs0[1] + (xs1[1]-xs0[1])*t;
+ c[0] = xb0[0] + (xb1[0]-xb0[0])*t;
+ c[1] = xb0[1] + (xb1[1]-xb0[1])*t;
+ double dtheta = theta1 - theta0;
+ double theta = theta0 + dtheta*t;
+ double cosT = cos(theta);
+ double sinT = sin(theta);
+
+ f = (p[1]-c[1]) * cosT - (p[0]-c[0]) * sinT;
+ df = ((xs1[1]-xs0[1]) - (xb1[1]-xb0[1]))*cosT - (p[1]-c[1])*sinT*dtheta -
+ ((xs1[0]-xs0[0]) - (xb1[0]-xb0[0]))*sinT - (p[0]-c[0])*cosT*dtheta;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixSRD::triside(double t, double &f, double &df)
+{
+ double p[2],c[2];
+
+ p[0] = xs0[0] + (xs1[0]-xs0[0])*t;
+ p[1] = xs0[1] + (xs1[1]-xs0[1])*t;
+ c[0] = xb0[0] + (xb1[0]-xb0[0])*t;
+ c[1] = xb0[1] + (xb1[1]-xb0[1])*t;
+ double dtheta = theta1 - theta0;
+ double theta = theta0 + dtheta*t;
+ double cosT = cos(theta);
+ double sinT = sin(theta);
+
+ f = (p[1]-c[1]) * cosT - (p[0]-c[0]) * sinT;
+ df = ((xs1[1]-xs0[1]) - (xb1[1]-xb0[1]))*cosT - (p[1]-c[1])*sinT*dtheta -
+ ((xs1[0]-xs0[0]) - (xb1[0]-xb0[0]))*sinT - (p[0]-c[0])*cosT*dtheta;
+}
+
+/* ---------------------------------------------------------------------- */
+
double FixSRD::memory_usage()
{
double bytes = 0.0;
diff --git a/src/SRD/fix_srd.h b/src/SRD/fix_srd.h
index 8552078e6d..43d9f7c2e3 100644
--- a/src/SRD/fix_srd.h
+++ b/src/SRD/fix_srd.h
@@ -43,9 +43,9 @@ class FixSRD : public Fix {
int me,nprocs;
int bigexist,biggroup,biggroupbit;
int collidestyle,lamdaflag,overlap,insideflag,exactflag,maxbounceallow;
- int cubicflag,shiftuser,shiftseed,shiftflag,streamflag;
+ int cubicflag,shiftuser,shiftseed,shiftflag,tstat;
double gridsrd,gridsearch,lamda,radfactor,cubictol;
- int triclinic,change_size,change_shape;
+ int triclinic,change_size,change_shape,deformflag;
double dt_big,dt_srd;
double mass_big,mass_srd;
@@ -64,6 +64,8 @@ class FixSRD : public Fix {
double walltrigger;
class AtomVecEllipsoid *avec_ellipsoid;
+ class AtomVecLine *avec_line;
+ class AtomVecTri *avec_tri;
// for orthogonal box, these are in box units
// for triclinic box, these are in lamda units
@@ -92,18 +94,21 @@ class FixSRD : public Fix {
struct Big {
int index; // local index of particle/wall
- int type; // SPHERE or ELLIPSOID or WALL
+ int type; // SPHERE or ELLIPSOID or LINE or TRI or WALL
double radius,radsq; // radius of sphere
double aradsqinv; // 3 ellipsoid radii
double bradsqinv;
double cradsqinv;
+ double length; // length of line segment
+ double normbody[3]; // normal of tri in body-frame
double cutbinsq; // add big to bin if within this distance
- double omega[3]; // current omega for sphere or ellipsoid
- double ex[3],ey[3],ez[3]; // current orientation vecs for ellipsoid
+ double omega[3]; // current omega for sphere/ellipsoid/tri/line
+ double ex[3],ey[3],ez[3]; // current orientation vecs for ellipsoid/tri
+ double norm[3]; // current unit normal of tri in space-frame
+ double theta; // current orientation of line
};
Big *biglist; // list of info for each owned & ghost big and wall
- int any_ellipsoids; // 1 if any big particles are ellipsoids
int torqueflag; // 1 if any big particle is torqued
// current size of particle-based arrays
@@ -123,7 +128,7 @@ class FixSRD : public Fix {
int owner; // 1 if I am owner of this bin, 0 if not
int n; // # of SRD particles in bin
double xctr[3]; // center point of bin, only used for triclinic
- double vave[3]; // sum of v components for SRD particles in bin
+ double vsum[3]; // sum of v components for SRD particles in bin
double random; // random value if I am owner
};
@@ -167,6 +172,17 @@ class FixSRD : public Fix {
int maxstencil; // max # of bins stencil array can hold
int **stencil; // list of 3d bin offsets a big particle can overlap
+ // persistent data for line/tri collision calculations
+
+ double tfraction,theta0,theta1;
+ double xs0[3],xs1[3],xsc[3];
+ double xb0[3],xb1[3],xbc[3];
+ double nbc[3];
+
+ // shared data for triangle collision calculations
+
+ // private functions
+
void reset_velocities();
void vbin_comm(int);
void vbin_pack(BinAve *, int, int *, double *);
@@ -177,6 +193,8 @@ class FixSRD : public Fix {
int inside_sphere(double *, double *, Big *);
int inside_ellipsoid(double *, double *, Big *);
+ int inside_line(double *, double *, double *, double *, Big *, double);
+ int inside_tri(double *, double *, double *, double *, Big *, double);
int inside_wall(double *, int);
double collision_sphere_exact(double *, double *, double *, double *,
@@ -187,18 +205,19 @@ class FixSRD : public Fix {
Big *, double *, double *, double *);
void collision_ellipsoid_inexact(double *, double *,
Big *, double *, double *, double *);
+ double collision_line_exact(double *, double *, double *, double *,
+ Big *, double, double *, double *, double *);
+ double collision_tri_exact(double *, double *, double *, double *,
+ Big *, double, double *, double *, double *);
double collision_wall_exact(double *, int, double *,
double *, double *, double *);
void collision_wall_inexact(double *, int, double *, double *, double *);
- void slip_sphere(double *, double *, double *, double *);
- void slip_ellipsoid(double *, double *, double *, Big *,
- double *, double *, double *);
+ void slip(double *, double *, double *, Big *,
+ double *, double *, double *);
void slip_wall(double *, int, double *, double *);
-
- void noslip(double *, double *, double *, Big *,
+ void noslip(double *, double *, double *, Big *, int,
double *, double *, double *);
- void noslip_wall(double *, int, double *, double *, double *);
void force_torque(double *, double *, double *,
double *, double *, double *);
@@ -217,11 +236,12 @@ class FixSRD : public Fix {
double point_bin_distance(double *, int, int, int);
double bin_bin_distance(int, int, int);
- void exyz_from_q(double *, double *, double *, double *);
- void omega_from_mq(double *, double *, double *, double *,
- double, double *, double *);
void velocity_stats(int);
+ double newton_raphson(double, double);
+ void lineside(double, double &, double &);
+ void triside(double, double &, double &);
+
double distance(int, int);
void print_collision(int, int, int, double, double,
double *, double *, double *, int);
diff --git a/src/atom.cpp b/src/atom.cpp
index f42b1250f1..023d7a173f 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -75,7 +75,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
radius = rmass = NULL;
vfrac = s0 = NULL;
x0 = NULL;
- ellipsoid = NULL;
+ ellipsoid = line = tri = NULL;
spin = NULL;
eradius = ervel = erforce = NULL;
cs = csforce = vforce = ervelforce = NULL;
@@ -106,7 +106,8 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
// initialize atom style and array existence flags
// customize by adding new flag
- sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
+ sphere_flag = ellipsoid_flag = line_flag = tri_flag = 0;
+ peri_flag = electron_flag = 0;
wavepacket_flag = sph_flag = 0;
molecule_flag = q_flag = mu_flag = 0;
@@ -185,6 +186,8 @@ Atom::~Atom()
memory->destroy(s0);
memory->destroy(x0);
memory->destroy(ellipsoid);
+ memory->destroy(line);
+ memory->destroy(tri);
memory->destroy(spin);
memory->destroy(eradius);
memory->destroy(ervel);
@@ -259,8 +262,9 @@ void Atom::create_avec(const char *style, int narg, char **arg, char *suffix)
// may have been set by old avec
// customize by adding new flag
- sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0;
-
+ sphere_flag = ellipsoid_flag = line_flag = tri_flag = 0;
+ peri_flag = electron_flag = 0;
+
molecule_flag = 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;
diff --git a/src/atom.h b/src/atom.h
index 98b687ec9a..429aacb1c2 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -51,7 +51,7 @@ class Atom : protected Pointers {
double **omega,**angmom,**torque;
double *radius,*rmass,*vfrac,*s0;
double **x0;
- int *ellipsoid;
+ int *ellipsoid,*line,*tri;
int *spin;
double *eradius,*ervel,*erforce,*ervelforce;
double *cs,*csforce,*vforce;
@@ -84,7 +84,7 @@ class Atom : protected Pointers {
// atom style and per-atom array existence flags
// customize by adding new flag
- int sphere_flag,ellipsoid_flag,peri_flag,electron_flag;
+ int sphere_flag,ellipsoid_flag,line_flag,tri_flag,peri_flag,electron_flag;
int wavepacket_flag,sph_flag;
int molecule_flag,q_flag,mu_flag;
diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp
index f7605d5552..80efd8a5bf 100644
--- a/src/compute_property_atom.cpp
+++ b/src/compute_property_atom.cpp
@@ -11,10 +11,14 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
+#include "math.h"
#include "string.h"
#include "compute_property_atom.h"
+#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
@@ -173,39 +177,39 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
pack_choice[i] = &ComputePropertyAtom::pack_angmomz;
} else if (strcmp(arg[iarg],"shapex") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapex;
} else if (strcmp(arg[iarg],"shapey") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapey;
} else if (strcmp(arg[iarg],"shapez") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_shapez;
} else if (strcmp(arg[iarg],"quatw") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatw;
} else if (strcmp(arg[iarg],"quati") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quati;
} else if (strcmp(arg[iarg],"quatj") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatj;
} else if (strcmp(arg[iarg],"quatk") == 0) {
- avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
- if (!avec) error->all(FLERR,"Compute property/atom for "
- "atom property that isn't allocated");
+ avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ if (!avec_ellipsoid) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_quatk;
} else if (strcmp(arg[iarg],"tqx") == 0) {
if (!atom->torque_flag)
@@ -244,6 +248,83 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
"atom property that isn't allocated");
pack_choice[i] = &ComputePropertyAtom::pack_erforce;
+ } else if (strcmp(arg[iarg],"end1x") == 0) {
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ if (!avec_line) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_end1x;
+ } else if (strcmp(arg[iarg],"end1y") == 0) {
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ if (!avec_line) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_end1y;
+ } else if (strcmp(arg[iarg],"end1z") == 0) {
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ if (!avec_line) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_end1z;
+ } else if (strcmp(arg[iarg],"end2x") == 0) {
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ if (!avec_line) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_end2x;
+ } else if (strcmp(arg[iarg],"end2y") == 0) {
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ if (!avec_line) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_end2y;
+ } else if (strcmp(arg[iarg],"end2z") == 0) {
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ if (!avec_line) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_end2z;
+
+ } else if (strcmp(arg[iarg],"corner1x") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner1x;
+ } else if (strcmp(arg[iarg],"corner1y") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner1y;
+ } else if (strcmp(arg[iarg],"corner1z") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner1z;
+ } else if (strcmp(arg[iarg],"corner2x") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner2x;
+ } else if (strcmp(arg[iarg],"corner2y") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner2y;
+ } else if (strcmp(arg[iarg],"corner2z") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner2z;
+ } else if (strcmp(arg[iarg],"corner3x") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner3x;
+ } else if (strcmp(arg[iarg],"corner3y") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner3y;
+ } else if (strcmp(arg[iarg],"corner3z") == 0) {
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
+ if (!avec_tri) error->all(FLERR,"Compute property/atom for "
+ "atom property that isn't allocated");
+ pack_choice[i] = &ComputePropertyAtom::pack_corner3z;
+
} else error->all(FLERR,"Invalid keyword in compute property/atom command");
}
@@ -994,7 +1075,7 @@ void ComputePropertyAtom::pack_angmomz(int n)
void ComputePropertyAtom::pack_shapex(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1011,7 +1092,7 @@ void ComputePropertyAtom::pack_shapex(int n)
void ComputePropertyAtom::pack_shapey(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1028,7 +1109,7 @@ void ComputePropertyAtom::pack_shapey(int n)
void ComputePropertyAtom::pack_shapez(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1045,7 +1126,7 @@ void ComputePropertyAtom::pack_shapez(int n)
void ComputePropertyAtom::pack_quatw(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1062,7 +1143,7 @@ void ComputePropertyAtom::pack_quatw(int n)
void ComputePropertyAtom::pack_quati(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1079,7 +1160,7 @@ void ComputePropertyAtom::pack_quati(int n)
void ComputePropertyAtom::pack_quatj(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1096,7 +1177,7 @@ void ComputePropertyAtom::pack_quatj(int n)
void ComputePropertyAtom::pack_quatk(int n)
{
- AtomVecEllipsoid::Bonus *bonus = avec->bonus;
+ AtomVecEllipsoid::Bonus *bonus = avec_ellipsoid->bonus;
int *ellipsoid = atom->ellipsoid;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@@ -1213,3 +1294,294 @@ void ComputePropertyAtom::pack_erforce(int n)
n += nvalues;
}
}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end1x(int n)
+{
+ AtomVecLine::Bonus *bonus = avec_line->bonus;
+ int *line = atom->line;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && line[i] >= 0)
+ buf[n] = x[i][0] - 0.5*bonus[line[i]].length*cos(bonus[line[i]].theta);
+ else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end1y(int n)
+{
+ AtomVecLine::Bonus *bonus = avec_line->bonus;
+ int *line = atom->line;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && line[i] >= 0)
+ buf[n] = x[i][1] - 0.5*bonus[line[i]].length*sin(bonus[line[i]].theta);
+ else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end1z(int n)
+{
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) buf[n] = x[i][2];
+ else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end2x(int n)
+{
+ AtomVecLine::Bonus *bonus = avec_line->bonus;
+ int *line = atom->line;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && line[i] >= 0)
+ buf[n] = x[i][0] + 0.5*bonus[line[i]].length*cos(bonus[line[i]].theta);
+ else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end2y(int n)
+{
+ AtomVecLine::Bonus *bonus = avec_line->bonus;
+ int *line = atom->line;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && line[i] >= 0)
+ buf[n] = x[i][1] + 0.5*bonus[line[i]].length*sin(bonus[line[i]].theta);
+ else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_end2z(int n)
+{
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) buf[n] = x[i][2];
+ else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner1x(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c1,c);
+ buf[n] = x[i][0] + c[0];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner1y(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c1,c);
+ buf[n] = x[i][1] + c[1];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner1z(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c1,c);
+ buf[n] = x[i][2] + c[2];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner2x(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c2,c);
+ buf[n] = x[i][0] + c[0];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner2y(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c2,c);
+ buf[n] = x[i][1] + c[1];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner2z(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c2,c);
+ buf[n] = x[i][2] + c[2];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner3x(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c3,c);
+ buf[n] = x[i][0] + c[0];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner3y(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c3,c);
+ buf[n] = x[i][1] + c[1];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputePropertyAtom::pack_corner3z(int n)
+{
+ AtomVecTri::Bonus *bonus = avec_tri->bonus;
+ int *tri = atom->tri;
+ double **x = atom->x;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ double p[3][3],c[3];
+ for (int i = 0; i < nlocal; i++) {
+ if ((mask[i] & groupbit) && tri[i] >= 0) {
+ MathExtra::quat_to_mat(bonus[tri[i]].quat,p);
+ MathExtra::matvec(p,bonus[tri[i]].c3,c);
+ buf[n] = x[i][2] + c[2];
+ } else buf[n] = 0.0;
+ n += nvalues;
+ }
+}
diff --git a/src/compute_property_atom.h b/src/compute_property_atom.h
index b32d0a93df..8f675b6204 100644
--- a/src/compute_property_atom.h
+++ b/src/compute_property_atom.h
@@ -38,7 +38,9 @@ class ComputePropertyAtom : public Compute {
double *vector;
double **array;
double *buf;
- class AtomVecEllipsoid *avec;
+ class AtomVecEllipsoid *avec_ellipsoid;
+ class AtomVecLine *avec_line;
+ class AtomVecTri *avec_tri;
typedef void (ComputePropertyAtom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
@@ -101,6 +103,21 @@ class ComputePropertyAtom : public Compute {
void pack_eradius(int);
void pack_ervel(int);
void pack_erforce(int);
+ void pack_end1x(int);
+ void pack_end1y(int);
+ void pack_end1z(int);
+ void pack_end2x(int);
+ void pack_end2y(int);
+ void pack_end2z(int);
+ void pack_corner1x(int);
+ void pack_corner1y(int);
+ void pack_corner1z(int);
+ void pack_corner2x(int);
+ void pack_corner2y(int);
+ void pack_corner2z(int);
+ void pack_corner3x(int);
+ void pack_corner3y(int);
+ void pack_corner3z(int);
};
}
diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index cf63dfa116..63df87b782 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -19,6 +19,8 @@
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
@@ -28,14 +30,20 @@
#include "random_mars.h"
#include "force.h"
#include "output.h"
+#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
+using namespace MathConst;
#define TOLERANCE 1.0e-6
#define EPSILON 1.0e-7
+#define SINERTIA 0.4 // moment of inertia prefactor for sphere
+#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid
+#define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment
+
/* ---------------------------------------------------------------------- */
FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
@@ -56,12 +64,12 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
// perform initial allocation of atom-based arrays
// register with Atom class
- extended = dorientflag = qorientflag = 0;
+ extended = orientflag = dorientflag = 0;
body = NULL;
displace = NULL;
eflags = NULL;
+ orient = NULL;
dorient = NULL;
- qorient = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
@@ -278,7 +286,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
else error->all(FLERR,"Illegal fix rigid command");
if (domain->dimension == 2 && (xflag == 1.0 || yflag == 1.0))
- error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation");
+ error->all(FLERR,"Fix rigid xy torque cannot be on for 2d simulation");
int count = 0;
for (int m = mlo; m <= mhi; m++) {
@@ -384,18 +392,24 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
// bitmasks for properties of extended particles
- INERTIA_POINT = 1;
- INERTIA_SPHERE = 2;
- INERTIA_ELLIPSOID = 4;
- ORIENT_DIPOLE = 8;
- ORIENT_QUAT = 16;
- OMEGA = 32;
- ANGMOM = 64;
- TORQUE = 128;
+ POINT = 1;
+ SPHERE = 2;
+ ELLIPSOID = 4;
+ LINE = 8;
+ TRIANGLE = 16;
+ DIPOLE = 32;
+ OMEGA = 64;
+ ANGMOM = 128;
+ TORQUE = 256;
+
+ MINUSPI = -MY_PI;
+ TWOPI = 2.0*MY_PI;
// atom style pointers to particles that store extra info
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
// print statistics
@@ -423,8 +437,8 @@ FixRigid::~FixRigid()
memory->destroy(body);
memory->destroy(displace);
memory->destroy(eflags);
+ memory->destroy(orient);
memory->destroy(dorient);
- memory->destroy(qorient);
// delete nbody-length arrays
@@ -469,7 +483,7 @@ int FixRigid::setmask()
void FixRigid::init()
{
- int i,ibody;
+ int i,itype,ibody;
triclinic = domain->triclinic;
@@ -504,24 +518,33 @@ void FixRigid::init()
// extended = 1 if any particle in a rigid body is finite size
// or has a dipole moment
- extended = dorientflag = qorientflag = 0;
+ extended = orientflag = dorientflag = 0;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
double **mu = atom->mu;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
int *type = atom->type;
int nlocal = atom->nlocal;
- if (atom->radius_flag || atom->ellipsoid_flag || atom->mu_flag) {
+ if (atom->radius_flag || atom->ellipsoid_flag || atom->line_flag ||
+ atom->tri_flag || atom->mu_flag) {
int flag = 0;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
if (radius && radius[i] > 0.0) flag = 1;
if (ellipsoid && ellipsoid[i] >= 0) flag = 1;
+ if (line && line[i] >= 0) flag = 1;
+ if (tri && tri[i] >= 0) flag = 1;
if (mu && mu[i][3] > 0.0) flag = 1;
}
@@ -529,35 +552,45 @@ void FixRigid::init()
}
// grow extended arrays and set extended flags for each particle
- // qorientflag = 1 if any particle stores quat orientation
+ // orientflag = 4 if any particle stores ellipsoid or tri orientation
+ // orientflag = 1 if any particle stores line orientation
// dorientflag = 1 if any particle stores dipole orientation
if (extended) {
+ if (atom->ellipsoid_flag) orientflag = 4;
+ if (atom->line_flag) orientflag = 1;
+ if (atom->tri_flag) orientflag = 4;
if (atom->mu_flag) dorientflag = 1;
- if (atom->ellipsoid_flag) qorientflag = 1;
grow_arrays(atom->nmax);
for (i = 0; i < nlocal; i++) {
eflags[i] = 0;
if (body[i] < 0) continue;
- // set INERTIA to POINT or SPHERE or ELLIPSOID
+ // set to POINT or SPHERE or ELLIPSOID or LINE
if (radius && radius[i] > 0.0) {
- eflags[i] |= INERTIA_SPHERE;
+ eflags[i] |= SPHERE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
} else if (ellipsoid && ellipsoid[i] >= 0) {
- eflags[i] |= INERTIA_ELLIPSOID;
- eflags[i] |= ORIENT_QUAT;
+ eflags[i] |= ELLIPSOID;
eflags[i] |= ANGMOM;
eflags[i] |= TORQUE;
- } else eflags[i] |= INERTIA_POINT;
+ } else if (line && line[i] >= 0) {
+ eflags[i] |= LINE;
+ eflags[i] |= OMEGA;
+ eflags[i] |= TORQUE;
+ } else if (tri && tri[i] >= 0) {
+ eflags[i] |= TRIANGLE;
+ eflags[i] |= ANGMOM;
+ eflags[i] |= TORQUE;
+ } else eflags[i] |= POINT;
// set DIPOLE if atom->mu and mu[3] > 0.0
if (atom->mu_flag && mu[i][3] > 0.0)
- eflags[i] |= ORIENT_DIPOLE;
+ eflags[i] |= DIPOLE;
}
}
@@ -592,8 +625,8 @@ void FixRigid::init()
if ((xbox && !periodicity[0]) || (ybox && !periodicity[1]) ||
(zbox && !periodicity[2]))
- error->one(FLERR,"Fix rigid atom has non-zero image flag "
- "in a non-periodic dimension");
+ error->one(FLERR,"Fix rigid atom has non-zero image flag "
+ "in a non-periodic dimension");
if (triclinic == 0) {
xunwrap = x[i][0] + xbox*xprd;
@@ -629,7 +662,7 @@ void FixRigid::init()
// dx,dy,dz = coords relative to center-of-mass
// symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
- double dx,dy,dz;
+ double dx,dy,dz,rad;
for (ibody = 0; ibody < nbody; ibody++)
for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
@@ -671,18 +704,19 @@ void FixRigid::init()
if (extended) {
double ivec[6];
- double *shape,*quatatom;
+ double *shape,*quatatom,*inertiaatom;
+ double length,theta;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
massone = rmass[i];
- if (eflags[i] & INERTIA_SPHERE) {
- sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
- sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
- sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
- } else if (eflags[i] & INERTIA_ELLIPSOID) {
+ if (eflags[i] & SPHERE) {
+ sum[ibody][0] += SINERTIA*massone * radius[i]*radius[i];
+ sum[ibody][1] += SINERTIA*massone * radius[i]*radius[i];
+ sum[ibody][2] += SINERTIA*massone * radius[i]*radius[i];
+ } else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
@@ -692,6 +726,26 @@ void FixRigid::init()
sum[ibody][3] += ivec[3];
sum[ibody][4] += ivec[4];
sum[ibody][5] += ivec[5];
+ } else if (eflags[i] & LINE) {
+ length = lbonus[line[i]].length;
+ theta = lbonus[line[i]].theta;
+ MathExtra::inertia_line(length,theta,massone,ivec);
+ sum[ibody][0] += ivec[0];
+ sum[ibody][1] += ivec[1];
+ sum[ibody][2] += ivec[2];
+ sum[ibody][3] += ivec[3];
+ sum[ibody][4] += ivec[4];
+ sum[ibody][5] += ivec[5];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::inertia_triangle(inertiaatom,quatatom,massone,ivec);
+ sum[ibody][0] += ivec[0];
+ sum[ibody][1] += ivec[1];
+ sum[ibody][2] += ivec[2];
+ sum[ibody][3] += ivec[3];
+ sum[ibody][4] += ivec[4];
+ sum[ibody][5] += ivec[5];
}
}
}
@@ -715,7 +769,8 @@ void FixRigid::init()
tensor[0][1] = tensor[1][0] = all[ibody][5];
ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
- if (ierror) error->all(FLERR,"Insufficient Jacobi rotations for rigid body");
+ if (ierror) error->all(FLERR,
+ "Insufficient Jacobi rotations for rigid body");
ex_space[ibody][0] = evectors[0][0];
ex_space[ibody][1] = evectors[1][0];
@@ -756,6 +811,7 @@ void FixRigid::init()
double qc[4],delta[3];
double *quatatom;
+ double theta_body;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) {
@@ -786,20 +842,34 @@ void FixRigid::init()
ez_space[ibody],delta,displace[i]);
if (extended) {
- if (eflags[i] & ORIENT_DIPOLE) {
+ if (eflags[i] & ELLIPSOID) {
+ quatatom = ebonus[ellipsoid[i]].quat;
+ MathExtra::qconjugate(quat[ibody],qc);
+ MathExtra::quatquat(qc,quatatom,orient[i]);
+ MathExtra::qnormalize(orient[i]);
+ } else if (eflags[i] & LINE) {
+ if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
+ else theta_body = -2.0*acos(quat[ibody][0]);
+ orient[i][0] = lbonus[line[i]].theta - theta_body;
+ while (orient[i][0] <= MINUSPI) orient[i][0] += TWOPI;
+ while (orient[i][0] > MY_PI) orient[i][0] -= TWOPI;
+ if (orientflag == 4) orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
+ } else if (eflags[i] & TRIANGLE) {
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::qconjugate(quat[ibody],qc);
+ MathExtra::quatquat(qc,quatatom,orient[i]);
+ MathExtra::qnormalize(orient[i]);
+ } else if (orientflag == 4) {
+ orient[i][0] = orient[i][1] = orient[i][2] = orient[i][3] = 0.0;
+ } else if (orientflag == 1)
+ orient[i][0] = 0.0;
+
+ if (eflags[i] & DIPOLE) {
MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],
ez_space[ibody],mu[i],dorient[i]);
MathExtra::snormalize3(mu[i][3],dorient[i],dorient[i]);
} else if (dorientflag)
dorient[i][0] = dorient[i][1] = dorient[i][2] = 0.0;
-
- if (eflags[i] & ORIENT_QUAT) {
- quatatom = ebonus[ellipsoid[i]].quat;
- MathExtra::qconjugate(quat[ibody],qc);
- MathExtra::quatquat(qc,quatatom,qorient[i]);
- MathExtra::qnormalize(qorient[i]);
- } else if (qorientflag)
- qorient[i][0] = qorient[i][1] = qorient[i][2] = qorient[i][3] = 0.0;
}
}
@@ -831,20 +901,39 @@ void FixRigid::init()
if (extended) {
double ivec[6];
- double *shape;
+ double *shape,*inertiaatom;
+ double length;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
massone = rmass[i];
- if (eflags[i] & INERTIA_SPHERE) {
- sum[ibody][0] += 0.4 * massone * radius[i]*radius[i];
- sum[ibody][1] += 0.4 * massone * radius[i]*radius[i];
- sum[ibody][2] += 0.4 * massone * radius[i]*radius[i];
- } else if (eflags[i] & INERTIA_ELLIPSOID) {
+ if (eflags[i] & SPHERE) {
+ sum[ibody][0] += SINERTIA*massone * radius[i]*radius[i];
+ sum[ibody][1] += SINERTIA*massone * radius[i]*radius[i];
+ sum[ibody][2] += SINERTIA*massone * radius[i]*radius[i];
+ } else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
- MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
+ MathExtra::inertia_ellipsoid(shape,orient[i],massone,ivec);
+ sum[ibody][0] += ivec[0];
+ sum[ibody][1] += ivec[1];
+ sum[ibody][2] += ivec[2];
+ sum[ibody][3] += ivec[3];
+ sum[ibody][4] += ivec[4];
+ sum[ibody][5] += ivec[5];
+ } else if (eflags[i] & LINE) {
+ length = lbonus[line[i]].length;
+ MathExtra::inertia_line(length,orient[i][0],massone,ivec);
+ sum[ibody][0] += ivec[0];
+ sum[ibody][1] += ivec[1];
+ sum[ibody][2] += ivec[2];
+ sum[ibody][3] += ivec[3];
+ sum[ibody][4] += ivec[4];
+ sum[ibody][5] += ivec[5];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ MathExtra::inertia_triangle(inertiaatom,orient[i],massone,ivec);
sum[ibody][0] += ivec[0];
sum[ibody][1] += ivec[1];
sum[ibody][2] += ivec[2];
@@ -894,7 +983,8 @@ void FixRigid::init()
ndof += fflag[ibody][0] + fflag[ibody][1] + fflag[ibody][2];
ndof += tflag[ibody][0] + tflag[ibody][1] + tflag[ibody][2];
}
- tfactor = force->mvv2e / (ndof * force->boltz);
+ if (ndof > 0.0) tfactor = force->mvv2e / (ndof * force->boltz);
+ else tfactor = 0.0;
}
/* ---------------------------------------------------------------------- */
@@ -996,24 +1086,29 @@ void FixRigid::setup(int vflag)
// extended particles add their rotation/torque to angmom/torque of body
if (extended) {
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **torque_one = atom->torque;
double *radius = atom->radius;
+ int *line = atom->line;
for (i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
if (eflags[i] & OMEGA) {
- if (rmass) massone = rmass[i];
- else massone = mass[type[i]];
- radone = radius[i];
- sum[ibody][0] += 0.4 * massone * radone*radone * omega_one[i][0];
- sum[ibody][1] += 0.4 * massone * radone*radone * omega_one[i][1];
- sum[ibody][2] += 0.4 * massone * radone*radone * omega_one[i][2];
+ if (eflags[i] & SPHERE) {
+ radone = radius[i];
+ sum[ibody][0] += SINERTIA*rmass[i] * radone*radone * omega_one[i][0];
+ sum[ibody][1] += SINERTIA*rmass[i] * radone*radone * omega_one[i][1];
+ sum[ibody][2] += SINERTIA*rmass[i] * radone*radone * omega_one[i][2];
+ } else if (eflags[i] & LINE) {
+ radone = lbonus[line[i]].length;
+ sum[ibody][2] += LINERTIA*rmass[i] * radone*radone * omega_one[i][2];
+ }
}
-
if (eflags[i] & ANGMOM) {
sum[ibody][0] += angmom_one[i][0];
sum[ibody][1] += angmom_one[i][1];
@@ -1516,7 +1611,7 @@ void FixRigid::deform(int flag)
void FixRigid::set_xv()
{
- int ibody;
+ int ibody,itype;
int xbox,ybox,zbox;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double xy,xz,yz;
@@ -1622,43 +1717,64 @@ void FixRigid::set_xv()
// set orientation, omega, angmom of each extended particle
if (extended) {
- double *shape,*quatatom;
+ double theta_body,theta;
+ double *shape,*quatatom,*inertiaatom;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecLine::Bonus *lbonus;
+ if (avec_line) lbonus = avec_line->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
double **mu = atom->mu;
int *ellipsoid = atom->ellipsoid;
+ int *line = atom->line;
+ int *tri = atom->tri;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
-
- if (eflags[i] & ORIENT_DIPOLE) {
- MathExtra::quat_to_mat(quat[ibody],p);
- MathExtra::matvec(p,dorient[i],mu[i]);
- MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
- }
- if (eflags[i] & ORIENT_QUAT) {
- quatatom = ebonus[ellipsoid[i]].quat;
- MathExtra::quatquat(quat[ibody],qorient[i],quatatom);
- MathExtra::qnormalize(quatatom);
- }
- if (eflags[i] & OMEGA) {
+
+ if (eflags[i] & SPHERE) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
- }
- if (eflags[i] & ANGMOM) {
+ } else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
- ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
- ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
- ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
+ MathExtra::quatquat(quat[ibody],orient[i],quatatom);
+ MathExtra::qnormalize(quatatom);
+ ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
+ ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
+ ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
angmom_one[i]);
+ } else if (eflags[i] & LINE) {
+ if (quat[ibody][3] >= 0.0) theta_body = 2.0*acos(quat[ibody][0]);
+ else theta_body = -2.0*acos(quat[ibody][0]);
+ theta = orient[i][0] + theta_body;
+ while (theta <= MINUSPI) theta += TWOPI;
+ while (theta > MY_PI) theta -= TWOPI;
+ lbonus[line[i]].theta = theta;
+ omega_one[i][0] = omega[ibody][0];
+ omega_one[i][1] = omega[ibody][1];
+ omega_one[i][2] = omega[ibody][2];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::quatquat(quat[ibody],orient[i],quatatom);
+ MathExtra::qnormalize(quatatom);
+ MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+ MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,
+ inertiaatom,angmom_one[i]);
+ }
+ if (eflags[i] & DIPOLE) {
+ MathExtra::quat_to_mat(quat[ibody],p);
+ MathExtra::matvec(p,dorient[i],mu[i]);
+ MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
}
}
}
@@ -1672,8 +1788,9 @@ void FixRigid::set_xv()
void FixRigid::set_v()
{
- int ibody;
+ int ibody,itype;
int xbox,ybox,zbox;
+ double dx,dy,dz;
double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
double xy,xz,yz;
double ione[3],exone[3],eyone[3],ezone[3],delta[3],vr[6];
@@ -1761,32 +1878,44 @@ void FixRigid::set_v()
// set omega, angmom of each extended particle
if (extended) {
- double *shape,*quatatom;
+ double *shape,*quatatom,*inertiaatom;
AtomVecEllipsoid::Bonus *ebonus;
if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
+ AtomVecTri::Bonus *tbonus;
+ if (avec_tri) tbonus = avec_tri->bonus;
double **omega_one = atom->omega;
double **angmom_one = atom->angmom;
int *ellipsoid = atom->ellipsoid;
+ int *tri = atom->tri;
for (int i = 0; i < nlocal; i++) {
if (body[i] < 0) continue;
ibody = body[i];
- if (eflags[i] & OMEGA) {
+ if (eflags[i] & SPHERE) {
omega_one[i][0] = omega[ibody][0];
omega_one[i][1] = omega[ibody][1];
omega_one[i][2] = omega[ibody][2];
- }
- if (eflags[i] & ANGMOM) {
+ } else if (eflags[i] & ELLIPSOID) {
shape = ebonus[ellipsoid[i]].shape;
quatatom = ebonus[ellipsoid[i]].quat;
- ione[0] = 0.2*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
- ione[1] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
- ione[2] = 0.2*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
+ ione[0] = EINERTIA*rmass[i] * (shape[1]*shape[1] + shape[2]*shape[2]);
+ ione[1] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[2]*shape[2]);
+ ione[2] = EINERTIA*rmass[i] * (shape[0]*shape[0] + shape[1]*shape[1]);
MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,ione,
angmom_one[i]);
+ } else if (eflags[i] & LINE) {
+ omega_one[i][0] = omega[ibody][0];
+ omega_one[i][1] = omega[ibody][1];
+ omega_one[i][2] = omega[ibody][2];
+ } else if (eflags[i] & TRIANGLE) {
+ inertiaatom = tbonus[tri[i]].inertia;
+ quatatom = tbonus[tri[i]].quat;
+ MathExtra::q_to_exyz(quatatom,exone,eyone,ezone);
+ MathExtra::omega_to_angmom(omega[ibody],exone,eyone,ezone,
+ inertiaatom,angmom_one[i]);
}
}
}
@@ -1804,8 +1933,8 @@ double FixRigid::memory_usage()
bytes += maxvatom*6 * sizeof(double);
if (extended) {
bytes += nmax * sizeof(int);
+ if (orientflag) bytes = nmax*orientflag * sizeof(double);
if (dorientflag) bytes = nmax*3 * sizeof(double);
- if (qorientflag) bytes = nmax*4 * sizeof(double);
}
return bytes;
}
@@ -1820,8 +1949,8 @@ void FixRigid::grow_arrays(int nmax)
memory->grow(displace,nmax,3,"rigid:displace");
if (extended) {
memory->grow(eflags,nmax,"rigid:eflags");
+ if (orientflag) memory->grow(orient,nmax,orientflag,"rigid:orient");
if (dorientflag) memory->grow(dorient,nmax,3,"rigid:dorient");
- if (qorientflag) memory->grow(qorient,nmax,4,"rigid:qorient");
}
}
@@ -1837,17 +1966,13 @@ void FixRigid::copy_arrays(int i, int j)
displace[j][2] = displace[i][2];
if (extended) {
eflags[j] = eflags[i];
+ for (int k = 0; k < orientflag; k++)
+ orient[j][k] = orient[i][k];
if (dorientflag) {
dorient[j][0] = dorient[i][0];
dorient[j][1] = dorient[i][1];
dorient[j][2] = dorient[i][2];
}
- if (qorientflag) {
- qorient[j][0] = qorient[i][0];
- qorient[j][1] = qorient[i][1];
- qorient[j][2] = qorient[i][2];
- qorient[j][3] = qorient[i][3];
- }
}
}
@@ -1877,17 +2002,13 @@ int FixRigid::pack_exchange(int i, double *buf)
int m = 4;
buf[m++] = eflags[i];
+ for (int j = 0; j < orientflag; j++)
+ buf[m++] = orient[i][j];
if (dorientflag) {
buf[m++] = dorient[i][0];
buf[m++] = dorient[i][1];
buf[m++] = dorient[i][2];
}
- if (qorientflag) {
- buf[m++] = qorient[i][0];
- buf[m++] = qorient[i][1];
- buf[m++] = qorient[i][2];
- buf[m++] = qorient[i][3];
- }
return m;
}
@@ -1905,17 +2026,13 @@ int FixRigid::unpack_exchange(int nlocal, double *buf)
int m = 4;
eflags[nlocal] = static_cast<int> (buf[m++]);
+ for (int j = 0; j < orientflag; j++)
+ orient[nlocal][j] = buf[m++];
if (dorientflag) {
dorient[nlocal][0] = buf[m++];
dorient[nlocal][1] = buf[m++];
dorient[nlocal][2] = buf[m++];
}
- if (qorientflag) {
- qorient[nlocal][0] = buf[m++];
- qorient[nlocal][1] = buf[m++];
- qorient[nlocal][2] = buf[m++];
- qorient[nlocal][3] = buf[m++];
- }
return m;
}
diff --git a/src/fix_rigid.h b/src/fix_rigid.h
index 06121ad47a..fa803df5f7 100644
--- a/src/fix_rigid.h
+++ b/src/fix_rigid.h
@@ -56,6 +56,7 @@ class FixRigid : public Fix {
double dtv,dtf,dtq;
double *step_respa;
int triclinic;
+ double MINUSPI,TWOPI;
int nbody; // # of rigid bodies
int *nrigid; // # of atoms in each rigid body
@@ -82,11 +83,11 @@ class FixRigid : public Fix {
int **remapflag; // PBC remap flags for each rigid body
int extended; // 1 if any particles have extended attributes
+ int orientflag; // 1 if particles store spatial orientation
int dorientflag; // 1 if particles store dipole orientation
- int qorientflag; // 1 if particles store quat orientation
int *eflags; // flags for extended particles
- double **qorient; // rotation state of ext particle wrt rigid body
+ double **orient; // orientation vector of particle wrt rigid body
double **dorient; // orientation of dipole mu wrt rigid body
double tfactor; // scale factor on temperature of rigid bodies
@@ -104,10 +105,10 @@ class FixRigid : public Fix {
class RanMars *random;
class AtomVecEllipsoid *avec_ellipsoid;
+ class AtomVecLine *avec_line;
+ class AtomVecTri *avec_tri;
- // bitmasks for eflags
- int INERTIA_POINT,INERTIA_SPHERE,INERTIA_ELLIPSOID;
- int ORIENT_DIPOLE,ORIENT_QUAT;
+ int POINT,SPHERE,ELLIPSOID,LINE,TRIANGLE,DIPOLE; // bitmasks for eflags
int OMEGA,ANGMOM,TORQUE;
void no_squish_rotate(int, double *, double *, double *, double);
diff --git a/src/read_data.cpp b/src/read_data.cpp
index a5a5e47ccb..772c3a3824 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -21,6 +21,8 @@
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
#include "comm.h"
#include "update.h"
#include "force.h"
@@ -42,7 +44,10 @@ using namespace LAMMPS_NS;
#define DELTA 4 // must be 2 or larger
// customize for new sections
-#define NSECTIONS 21 // change when add to header::section_keywords
+#define NSECTIONS 23 // change when add to header::section_keywords
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
/* ---------------------------------------------------------------------- */
@@ -60,6 +65,10 @@ ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp)
nellipsoids = 0;
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ nlines = 0;
+ avec_line = (AtomVecLine *) atom->style_match("line");
+ ntris = 0;
+ avec_tri = (AtomVecTri *) atom->style_match("tri");
}
/* ---------------------------------------------------------------------- */
@@ -151,7 +160,17 @@ void ReadData::command(int narg, char **arg)
if (!avec_ellipsoid)
error->all(FLERR,"Invalid data file section: Ellipsoids");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Ellipsoids");
- ellipsoids();
+ bonus(nellipsoids,(AtomVec *) avec_ellipsoid,"ellipsoids");
+ } else if (strcmp(keyword,"Lines") == 0) {
+ if (!avec_line)
+ error->all(FLERR,"Invalid data file section: Lines");
+ if (atomflag == 0) error->all(FLERR,"Must read Atoms before Lines");
+ bonus(nlines,(AtomVec *) avec_line,"lines");
+ } else if (strcmp(keyword,"Triangles") == 0) {
+ if (!avec_tri)
+ error->all(FLERR,"Invalid data file section: Triangles");
+ if (atomflag == 0) error->all(FLERR,"Must read Atoms before Triangles");
+ bonus(ntris,(AtomVec *) avec_tri,"triangles");
} else if (strcmp(keyword,"Bonds") == 0) {
if (atom->avec->bonds_allow == 0)
@@ -222,25 +241,31 @@ void ReadData::command(int narg, char **arg)
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before MiddleBondTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before "
+ "MiddleBondTorsion Coeffs");
dihedralcoeffs(1);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before EndBondTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before EndBondTorsion Coeffs");
dihedralcoeffs(2);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before AngleTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before AngleTorsion Coeffs");
dihedralcoeffs(3);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before AngleAngleTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before "
+ "AngleAngleTorsion Coeffs");
dihedralcoeffs(4);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
@@ -274,7 +299,8 @@ void ReadData::command(int narg, char **arg)
// error if natoms > 0 yet no atoms were read
- if (atom->natoms > 0 && atomflag == 0) error->all(FLERR,"No atoms in data file");
+ if (atom->natoms > 0 && atomflag == 0)
+ error->all(FLERR,"No atoms in data file");
// create bond topology now that system is defined
@@ -303,7 +329,7 @@ void ReadData::header(int flag)
// customize for new sections
char *section_keywords[NSECTIONS] =
- {"Atoms","Velocities","Ellipsoids",
+ {"Atoms","Velocities","Ellipsoids","Lines","Triangles",
"Bonds","Angles","Dihedrals","Impropers",
"Masses","Pair Coeffs","Bond Coeffs","Angle Coeffs",
"Dihedral Coeffs","Improper Coeffs",
@@ -373,6 +399,14 @@ void ReadData::header(int flag)
if (!avec_ellipsoid)
error->all(FLERR,"No ellipsoids allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nellipsoids);
+ } else if (strstr(line,"lines")) {
+ if (!avec_line)
+ error->all(FLERR,"No lines allowed with this atom style");
+ sscanf(line,BIGINT_FORMAT,&nlines);
+ } else if (strstr(line,"triangles")) {
+ if (!avec_tri)
+ error->all(FLERR,"No triangles allowed with this atom style");
+ sscanf(line,BIGINT_FORMAT,&ntris);
}
else if (strstr(line,"xlo xhi"))
@@ -475,7 +509,8 @@ void ReadData::atoms()
if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms);
}
- if (natoms != atom->natoms) error->all(FLERR,"Did not assign all atoms correctly");
+ if (natoms != atom->natoms)
+ error->all(FLERR,"Did not assign all atoms correctly");
// if any atom ID < 0, error
// if all atom IDs = 0, tag_enable = 0
@@ -568,11 +603,11 @@ void ReadData::velocities()
}
/* ----------------------------------------------------------------------
- read all ellipsoids
+ read all bonus data
to find atoms, must build atom map if not a molecular system
------------------------------------------------------------------------- */
-void ReadData::ellipsoids()
+void ReadData::bonus(bigint nbonus, AtomVec *ptr, char *type)
{
int i,m,nchunk;
@@ -585,7 +620,7 @@ void ReadData::ellipsoids()
}
bigint nread = 0;
- bigint natoms = nellipsoids;
+ bigint natoms = nbonus;
while (nread < natoms) {
if (natoms-nread > CHUNK) nchunk = CHUNK;
@@ -603,7 +638,7 @@ void ReadData::ellipsoids()
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
- atom->data_bonus(nchunk,buffer,avec_ellipsoid);
+ atom->data_bonus(nchunk,buffer,ptr);
nread += nchunk;
}
@@ -613,8 +648,8 @@ void ReadData::ellipsoids()
}
if (me == 0) {
- if (screen) fprintf(screen," " BIGINT_FORMAT " ellipsoids\n",natoms);
- if (logfile) fprintf(logfile," " BIGINT_FORMAT " ellipsoids\n",natoms);
+ if (screen) fprintf(screen," " BIGINT_FORMAT " %s\n",natoms,type);
+ if (logfile) fprintf(logfile," " BIGINT_FORMAT " %s\n",natoms,type);
}
}
@@ -660,7 +695,8 @@ void ReadData::bonds()
if (screen) fprintf(screen," " BIGINT_FORMAT " bonds\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",sum/factor);
}
- if (sum != factor*atom->nbonds) error->all(FLERR,"Bonds assigned incorrectly");
+ if (sum != factor*atom->nbonds)
+ error->all(FLERR,"Bonds assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
@@ -705,7 +741,8 @@ void ReadData::angles()
if (screen) fprintf(screen," " BIGINT_FORMAT " angles\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n",sum/factor);
}
- if (sum != factor*atom->nangles) error->all(FLERR,"Angles assigned incorrectly");
+ if (sum != factor*atom->nangles)
+ error->all(FLERR,"Angles assigned incorrectly");
}
/* ---------------------------------------------------------------------- */
@@ -1010,6 +1047,8 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
int natoms = static_cast<int> (atom->natoms);
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
int ellipsoid_flag = 0;
+ int line_flag = 0;
+ int tri_flag = 0;
// customize for new sections
// allocate topology counting vector
@@ -1031,6 +1070,14 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
error->all(FLERR,"Invalid data file section: Ellipsoids");
ellipsoid_flag = 1;
skip_lines(nellipsoids);
+ } else if (strcmp(keyword,"Lines") == 0) {
+ if (!avec_line) error->all(FLERR,"Invalid data file section: Lines");
+ line_flag = 1;
+ skip_lines(nlines);
+ } else if (strcmp(keyword,"Triangles") == 0) {
+ if (!avec_tri) error->all(FLERR,"Invalid data file section: Triangles");
+ tri_flag = 1;
+ skip_lines(ntris);
} else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
@@ -1077,25 +1124,31 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before MiddleBondTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before "
+ "MiddleBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before EndBondTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before EndBondTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before AngleTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before AngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
- error->all(FLERR,"Must define dihedral_style before AngleAngleTorsion Coeffs");
+ error->all(FLERR,
+ "Must define dihedral_style before "
+ "AngleAngleTorsion Coeffs");
skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
@@ -1252,6 +1305,10 @@ void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
if (nellipsoids && !ellipsoid_flag)
error->one(FLERR,"Needed bonus data not in data file");
+ if (nlines && !line_flag)
+ error->one(FLERR,"Needed bonus data not in data file");
+ if (ntris && !tri_flag)
+ error->one(FLERR,"Needed bonus data not in data file");
}
/* ----------------------------------------------------------------------
diff --git a/src/read_data.h b/src/read_data.h
index 80fd441413..67398469c5 100644
--- a/src/read_data.h
+++ b/src/read_data.h
@@ -40,6 +40,10 @@ class ReadData : protected Pointers {
bigint nellipsoids;
class AtomVecEllipsoid *avec_ellipsoid;
+ bigint nlines;
+ class AtomVecLine *avec_line;
+ bigint ntris;
+ class AtomVecTri *avec_tri;
void open(char *);
void scan(int &, int &, int &, int &);
@@ -51,7 +55,7 @@ class ReadData : protected Pointers {
void atoms();
void velocities();
- void ellipsoids();
+ void bonus(bigint, class AtomVec *, char *);
void bonds();
void angles();
diff --git a/src/set.cpp b/src/set.cpp
index 4e76dfa2b4..2a4c25731b 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -19,6 +19,8 @@
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
+#include "atom_vec_line.h"
+#include "atom_vec_tri.h"
#include "domain.h"
#include "region.h"
#include "group.h"
@@ -35,8 +37,8 @@ using namespace LAMMPS_NS;
using namespace MathConst;
enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT};
-enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,
- DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,
+enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI,
+ DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM,THETA,ANGMOM,
DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER,
MESO_E,MESO_CV,MESO_RHO};
@@ -150,6 +152,22 @@ void Set::command(int narg, char **arg)
}
set(SHAPE);
iarg += 4;
+ } else if (strcmp(arg[iarg],"length") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+ dvalue = atof(arg[iarg+1]);
+ if (!atom->line_flag)
+ error->all(FLERR,"Cannot set this attribute for this atom style");
+ if (dvalue < 0.0) error->all(FLERR,"Invalid length in set command");
+ set(LENGTH);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"tri") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+ dvalue = atof(arg[iarg+1]);
+ if (!atom->tri_flag)
+ error->all(FLERR,"Cannot set this attribute for this atom style");
+ if (dvalue < 0.0) error->all(FLERR,"Invalid length in set command");
+ set(TRI);
+ iarg += 2;
} else if (strcmp(arg[iarg],"dipole") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal set command");
xvalue = atof(arg[iarg+1]);
@@ -177,19 +195,36 @@ void Set::command(int narg, char **arg)
yvalue = atof(arg[iarg+2]);
zvalue = atof(arg[iarg+3]);
wvalue = atof(arg[iarg+4]);
- if (!atom->ellipsoid_flag)
+ if (!atom->ellipsoid_flag && !atom->tri_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)
+ if (!atom->ellipsoid_flag && !atom->tri_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
if (ivalue <= 0)
error->all(FLERR,"Invalid random number seed in set command");
setrandom(QUAT_RANDOM);
iarg += 2;
+ } else if (strcmp(arg[iarg],"theta") == 0) {
+ if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
+ dvalue = atof(arg[iarg+1]);
+ dvalue *= MY_PI/180.0;
+ if (!atom->line_flag)
+ error->all(FLERR,"Cannot set this attribute for this atom style");
+ set(THETA);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"angmom") == 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 && !atom->tri_flag)
+ error->all(FLERR,"Cannot set this attribute for this atom style");
+ set(ANGMOM);
+ iarg += 4;
} else if (strcmp(arg[iarg],"diameter") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
dvalue = atof(arg[iarg+1]);
@@ -385,6 +420,8 @@ void Set::set(int keyword)
{
AtomVecEllipsoid *avec_ellipsoid =
(AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
+ AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
selection(atom->nlocal);
@@ -405,23 +442,49 @@ void Set::set(int keyword)
else if (keyword == MESO_CV) atom->cv[i] = dvalue;
else if (keyword == MESO_RHO) atom->rho[i] = dvalue;
- // set shape
+ // set shape of ellipsoidal particle
else if (keyword == SHAPE)
avec_ellipsoid->set_shape(i,0.5*xvalue,0.5*yvalue,0.5*zvalue);
+ // set length of line particle
+
+ else if (keyword == LENGTH)
+ avec_line->set_length(i,dvalue);
+
+ // set corners of tri particle
+
+ else if (keyword == TRI)
+ avec_tri->set_equilateral(i,dvalue);
+
// set rmass via density
// if radius > 0.0, treat as sphere
// if shape > 0.0, treat as ellipsoid
+ // if length > 0.0, treat as line
+ // if area > 0.0, treat as tri
// else set rmass to density directly
else if (keyword == DENSITY) {
if (atom->radius_flag && atom->radius[i] > 0.0)
- atom->rmass[i] = 4.0*MY_PI*THIRD *
+ atom->rmass[i] = 4.0*MY_PI/3.0 *
atom->radius[i]*atom->radius[i]*atom->radius[i] * dvalue;
else if (atom->ellipsoid_flag && atom->ellipsoid[i] >= 0) {
double *shape = avec_ellipsoid->bonus[atom->ellipsoid[i]].shape;
- atom->rmass[i] = 4.0*MY_PI*THIRD * shape[0]*shape[1]*shape[2] * dvalue;
+ atom->rmass[i] = 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2] * dvalue;
+ } else if (atom->line_flag && atom->line[i] >= 0) {
+ double length = avec_line->bonus[atom->line[i]].length;
+ atom->rmass[i] = length * dvalue;
+ } else if (atom->tri_flag && atom->tri[i] >= 0) {
+ double *c1 = avec_tri->bonus[atom->tri[i]].c1;
+ double *c2 = avec_tri->bonus[atom->tri[i]].c2;
+ double *c3 = avec_tri->bonus[atom->tri[i]].c3;
+ double c2mc1[2],c3mc1[3];
+ MathExtra::sub3(c2,c1,c2mc1);
+ MathExtra::sub3(c3,c1,c3mc1);
+ double norm[3];
+ MathExtra::cross3(c2mc1,c3mc1,norm);
+ double area = 0.5 * MathExtra::len3(norm);
+ atom->rmass[i] = area * dvalue;
} else atom->rmass[i] = dvalue;
// reset any or all of 3 image flags
@@ -446,13 +509,17 @@ void Set::set(int keyword)
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
+ // set quaternion orientation of ellipsoid or tri particle
} 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 *quat;
+ if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
+ quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
+ else if (avec_tri && atom->tri[i] >= 0)
+ quat = avec_tri->bonus[atom->tri[i]].quat;
+ else
+ error->one(FLERR,"Cannot set quaternion for atom that has none");
+
double theta2 = MY_PI2 * wvalue/180.0;
double sintheta2 = sin(theta2);
quat[0] = cos(theta2);
@@ -460,7 +527,22 @@ void Set::set(int keyword)
quat[2] = yvalue * sintheta2;
quat[3] = zvalue * sintheta2;
MathExtra::qnormalize(quat);
+
+ // set theta of line particle
+
+ } else if (keyword == THETA) {
+ if (atom->line[i] < 0)
+ error->one(FLERR,"Cannot set theta for atom that is not a line");
+ avec_line->bonus[atom->line[i]].theta = dvalue;
+
+ // set angmom of ellipsoidal or tri particle
+
+ } else if (keyword == ANGMOM) {
+ atom->angmom[i][0] = xvalue;
+ atom->angmom[i][1] = yvalue;
+ atom->angmom[i][2] = zvalue;
}
+
count++;
}
}
@@ -475,6 +557,11 @@ void Set::setrandom(int keyword)
{
int i;
+ AtomVecEllipsoid *avec_ellipsoid =
+ (AtomVecEllipsoid *) atom->style_match("ellipsoid");
+ AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
+ AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
+
selection(atom->nlocal);
RanPark *random = new RanPark(lmp,1);
double **x = atom->x;
@@ -535,13 +622,10 @@ void Set::setrandom(int keyword)
}
// 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 *tri = atom->tri;
int nlocal = atom->nlocal;
double *quat;
@@ -549,10 +633,13 @@ void Set::setrandom(int keyword)
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;
+ if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
+ quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
+ else if (avec_tri && atom->tri[i] >= 0)
+ quat = avec_tri->bonus[atom->tri[i]].quat;
+ else
+ error->one(FLERR,"Cannot set quaternion for atom that has none");
+
random->reset(seed,x[i]);
s = random->uniform();
t1 = sqrt(1.0-s);
@@ -570,10 +657,11 @@ void Set::setrandom(int keyword)
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;
+ if (avec_ellipsoid && atom->ellipsoid[i] >= 0)
+ quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat;
+ else
+ error->one(FLERR,"Cannot set quaternion for atom that has none");
+
random->reset(seed,x[i]);
theta2 = MY_PI*random->uniform();
quat[0] = cos(theta2);
--
GitLab