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

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

Browse Source
This commit is contained in:
sjplimp
2012-10-25 16:24:37 +00:00
parent cb9d45386c
commit 6f4bba7edf
11 changed files with 149 additions and 287 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/compute_com_molecule.cpp
View File

@ -96,8 +96,8 @@ void ComputeCOMMolecule::init()
void ComputeCOMMolecule::compute_array() void ComputeCOMMolecule::compute_array()
{ {
int i,imol; int i,imol;
double xbox,ybox,zbox;
double massone; double massone;
double unwrap[3];
invoked_array = update->ntimestep; invoked_array = update->ntimestep;
@ -113,23 +113,17 @@ void ComputeCOMMolecule::compute_array()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
imol = molecule[i]; imol = molecule[i];
if (molmap) imol = molmap[imol-idlo]; if (molmap) imol = molmap[imol-idlo];
else imol--; else imol--;
com[imol][0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
com[imol][1] += (x[i][1] + ybox*yprd) * massone; com[imol][0] += unwrap[0] * massone;
com[imol][2] += (x[i][2] + zbox*zprd) * massone; com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
} }
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,

17
src/compute_gyration.cpp
View File

@ -83,25 +83,22 @@ void ComputeGyration::compute_vector()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double rg[6]; double rg[6];
rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0; rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
rg[0] += dx*dx * massone; rg[0] += dx*dx * massone;
rg[1] += dy*dy * massone; rg[1] += dy*dy * massone;
rg[2] += dz*dz * massone; rg[2] += dz*dz * massone;

57
src/compute_gyration_molecule.cpp
View File

@ -123,8 +123,8 @@ void ComputeGyrationMolecule::init()
void ComputeGyrationMolecule::compute_vector() void ComputeGyrationMolecule::compute_vector()
{ {
int i,imol; int i,imol;
double xbox,ybox,zbox,dx,dy,dz; double dx,dy,dz,massone;
double massone; double unwrap[3];
invoked_array = update->ntimestep; invoked_array = update->ntimestep;
@ -141,21 +141,16 @@ void ComputeGyrationMolecule::compute_vector()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
imol = molecule[i]; imol = molecule[i];
if (molmap) imol = molmap[imol-idlo]; if (molmap) imol = molmap[imol-idlo];
else imol--; else imol--;
dx = (x[i][0] + xbox*xprd) - comall[imol][0]; domain->unmap(x[i],image[i],unwrap);
dy = (x[i][1] + ybox*yprd) - comall[imol][1]; dx = unwrap[0] - comall[imol][0];
dz = (x[i][2] + zbox*zprd) - comall[imol][2]; dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
rg[imol] += (dx*dx + dy*dy + dz*dz) * massone; rg[imol] += (dx*dx + dy*dy + dz*dz) * massone;
@ -171,8 +166,8 @@ void ComputeGyrationMolecule::compute_vector()
void ComputeGyrationMolecule::compute_array() void ComputeGyrationMolecule::compute_array()
{ {
int i,j,imol; int i,j,imol;
double xbox,ybox,zbox,dx,dy,dz; double dx,dy,dz,massone;
double massone; double unwrap[3];
invoked_array = update->ntimestep; invoked_array = update->ntimestep;
@ -191,21 +186,15 @@ void ComputeGyrationMolecule::compute_array()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
imol = molecule[i]; imol = molecule[i];
if (molmap) imol = molmap[imol-idlo]; if (molmap) imol = molmap[imol-idlo];
else imol--; else imol--;
dx = (x[i][0] + xbox*xprd) - comall[imol][0]; domain->unmap(x[i],image[i],unwrap);
dy = (x[i][1] + ybox*yprd) - comall[imol][1]; dx = unwrap[0] - comall[imol][0];
dz = (x[i][2] + zbox*zprd) - comall[imol][2]; dy = unwrap[1] - comall[imol][1];
dz = unwrap[2] - comall[imol][2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
rgt[imol][0] += dx*dx * massone; rgt[imol][0] += dx*dx * massone;
@ -233,8 +222,8 @@ void ComputeGyrationMolecule::compute_array()
void ComputeGyrationMolecule::molcom() void ComputeGyrationMolecule::molcom()
{ {
int i,imol; int i,imol;
double xbox,ybox,zbox,dx,dy,dz; double dx,dy,dz,massone;
double massone; double unwrap[3];
for (i = 0; i < nmolecules; i++) for (i = 0; i < nmolecules; i++)
com[i][0] = com[i][1] = com[i][2] = 0.0; com[i][0] = com[i][1] = com[i][2] = 0.0;
@ -248,23 +237,17 @@ void ComputeGyrationMolecule::molcom()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i]; imol = molecule[i];
if (molmap) imol = molmap[imol-idlo]; if (molmap) imol = molmap[imol-idlo];
else imol--; else imol--;
com[imol][0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
com[imol][1] += (x[i][1] + ybox*yprd) * massone; if (rmass) massone = rmass[i];
com[imol][2] += (x[i][2] + zbox*zprd) * massone; else massone = mass[type[i]];
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
} }
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,

22
src/compute_msd_molecule.cpp
View File

@ -111,8 +111,8 @@ void ComputeMSDMolecule::init()
void ComputeMSDMolecule::compute_array() void ComputeMSDMolecule::compute_array()
{ {
int i,imol; int i,imol;
double xbox,ybox,zbox,dx,dy,dz; double dx,dy,dz,massone;
double massone; double unwrap[3];
invoked_array = update->ntimestep; invoked_array = update->ntimestep;
@ -130,23 +130,17 @@ void ComputeMSDMolecule::compute_array()
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
imol = molecule[i]; imol = molecule[i];
if (molmap) imol = molmap[imol-idlo]; if (molmap) imol = molmap[imol-idlo];
else imol--; else imol--;
com[imol][0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
com[imol][1] += (x[i][1] + ybox*yprd) * massone; if (rmass) massone = rmass[i];
com[imol][2] += (x[i][2] + zbox*zprd) * massone; else massone = mass[type[i]];
com[imol][0] += unwrap[0] * massone;
com[imol][1] += unwrap[1] * massone;
com[imol][2] += unwrap[2] * massone;
} }
MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules, MPI_Allreduce(&com[0][0],&comall[0][0],3*nmolecules,

15
src/fix_addforce.cpp
View File

@ -217,12 +217,6 @@ void FixAddForce::min_setup(int vflag)
void FixAddForce::post_force(int vflag) void FixAddForce::post_force(int vflag)
{ {
int xbox,ybox,zbox;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
int *mask = atom->mask; int *mask = atom->mask;
@ -247,18 +241,15 @@ void FixAddForce::post_force(int vflag)
// potential energy = - x dot f in unwrapped coords // potential energy = - x dot f in unwrapped coords
if (varflag == CONSTANT) { if (varflag == CONSTANT) {
double unwrap[3];
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
if (iregion >= 0 && if (iregion >= 0 &&
!domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) !domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]))
continue; continue;
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; foriginal[0] -= xvalue*unwrap[0] + yvalue*unwrap[1] + zvalue*unwrap[2];
zbox = (image[i] >> IMG2BITS) - IMGMAX;
foriginal[0] -= xvalue * (x[i][0]+xbox*xprd) +
yvalue * (x[i][1]+ybox*yprd) + zvalue * (x[i][2]+zbox*zprd);
foriginal[1] += f[i][0]; foriginal[1] += f[i][0];
foriginal[2] += f[i][1]; foriginal[2] += f[i][1];
foriginal[3] += f[i][2]; foriginal[3] += f[i][2];

18
src/fix_momentum.cpp
View File

@ -57,7 +57,8 @@ FixMomentum::FixMomentum(LAMMPS *lmp, int narg, char **arg) :
if (linear) if (linear)
if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 || if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 ||
zflag < 0 || zflag > 1) error->all(FLERR,"Illegal fix momentum command"); zflag < 0 || zflag > 1)
error->all(FLERR,"Illegal fix momentum command");
// cannot have 0 atoms in group // cannot have 0 atoms in group
@ -121,20 +122,15 @@ void FixMomentum::end_of_step()
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - xcm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - xcm[1];
dx = (x[i][0] + xbox*xprd) - xcm[0]; dz = unwrap[2] - xcm[2];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
v[i][0] -= omega[1]*dz - omega[2]*dy; v[i][0] -= omega[1]*dz - omega[2]*dy;
v[i][1] -= omega[2]*dx - omega[0]*dz; v[i][1] -= omega[2]*dx - omega[0]*dz;
v[i][2] -= omega[0]*dy - omega[1]*dx; v[i][2] -= omega[0]*dy - omega[1]*dx;

15
src/fix_spring_rg.cpp
View File

@ -110,20 +110,15 @@ void FixSpringRG::post_force(int vflag)
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double massfrac; double massfrac;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
term1 = 2.0 * k * (1.0 - rg0/rg); term1 = 2.0 * k * (1.0 - rg0/rg);
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
dx = (x[i][0] + xbox*xprd) - xcm[0];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
massfrac = mass[type[i]]/masstotal; massfrac = mass[type[i]]/masstotal;
f[i][0] -= term1*dx*massfrac; f[i][0] -= term1*dx*massfrac;
f[i][1] -= term1*dy*massfrac; f[i][1] -= term1*dy*massfrac;

34
src/fix_spring_self.cpp
View File

@ -45,9 +45,10 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command"); if (k <= 0.0) error->all(FLERR,"Illegal fix spring/self command");
xflag = yflag = zflag = 1; xflag = yflag = zflag = 1;
if (narg == 5) { if (narg == 5) {
if (strcmp(arg[4],"xyz") == 0) { if (strcmp(arg[4],"xyz") == 0) {
; /* default */ xflag = yflag = zflag = 1;
} else if (strcmp(arg[4],"xy") == 0) { } else if (strcmp(arg[4],"xy") == 0) {
zflag = 0; zflag = 0;
} else if (strcmp(arg[4],"xz") == 0) { } else if (strcmp(arg[4],"xz") == 0) {
@ -78,20 +79,9 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) domain->unmap(x[i],image[i],xoriginal[i]);
xbox = (image[i] & IMGMASK) - IMGMAX; else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xoriginal[i][0] = x[i][0] + xbox*xprd;
xoriginal[i][1] = x[i][1] + ybox*yprd;
xoriginal[i][2] = x[i][2] + zbox*zprd;
} else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
} }
espring = 0.0; espring = 0.0;
@ -161,21 +151,17 @@ void FixSpringSelf::post_force(int vflag)
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double unwrap[3];
espring = 0.0; espring = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - xoriginal[i][0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - xoriginal[i][1];
dx = x[i][0] + xbox*xprd - xoriginal[i][0]; dz = unwrap[2] - xoriginal[i][2];
dy = x[i][1] + ybox*yprd - xoriginal[i][1];
dz = x[i][2] + zbox*zprd - xoriginal[i][2];
if (!xflag) dx = 0.0; if (!xflag) dx = 0.0;
if (!yflag) dy = 0.0; if (!yflag) dy = 0.0;
if (!zflag) dz = 0.0; if (!zflag) dz = 0.0;

51
src/fix_tmd.cpp
View File

@ -98,26 +98,17 @@ FixTMD::FixTMD(LAMMPS *lmp, int narg, char **arg) :
double *mass = atom->mass; double *mass = atom->mass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double dx,dy,dz; double dx,dy,dz;
int xbox,ybox,zbox;
rho_start = 0.0; rho_start = 0.0;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],xold[i]);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = xold[i][0] - xf[i][0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = xold[i][1] - xf[i][1];
dx = x[i][0] + xbox*xprd - xf[i][0]; dz = xold[i][2] - xf[i][2];
dy = x[i][1] + ybox*yprd - xf[i][1];
dz = x[i][2] + zbox*zprd - xf[i][2];
rho_start += mass[type[i]]*(dx*dx + dy*dy + dz*dz); rho_start += mass[type[i]]*(dx*dx + dy*dy + dz*dz);
xold[i][0] = x[i][0] + xbox*xprd;
xold[i][1] = x[i][1] + ybox*yprd;
xold[i][2] = x[i][2] + zbox*zprd;
} else xold[i][0] = xold[i][1] = xold[i][2] = 0.0; } else xold[i][0] = xold[i][1] = xold[i][2] = 0.0;
} }
@ -190,10 +181,8 @@ void FixTMD::initial_integrate(int vflag)
double dxold,dyold,dzold,xback,yback,zback; double dxold,dyold,dzold,xback,yback,zback;
double gamma_forward,gamma_back,gamma_max,lambda; double gamma_forward,gamma_back,gamma_max,lambda;
double kt,fr,kttotal,frtotal,dtfm; double kt,fr,kttotal,frtotal,dtfm;
double unwrap[3];
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -202,7 +191,6 @@ void FixTMD::initial_integrate(int vflag)
int *type = atom->type; int *type = atom->type;
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double delta = update->ntimestep - update->beginstep; double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep; delta /= update->endstep - update->beginstep;
@ -216,12 +204,10 @@ void FixTMD::initial_integrate(int vflag)
dxold = xold[i][0] - xf[i][0]; dxold = xold[i][0] - xf[i][0];
dyold = xold[i][1] - xf[i][1]; dyold = xold[i][1] - xf[i][1];
dzold = xold[i][2] - xf[i][2]; dzold = xold[i][2] - xf[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - xf[i][0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - xf[i][1];
dx = x[i][0] + xbox*xprd - xf[i][0]; dz = unwrap[2] - xf[i][2];
dy = x[i][1] + ybox*yprd - xf[i][1];
dz = x[i][2] + zbox*zprd - xf[i][2];
a += mass[type[i]]*(dxold*dxold + dyold*dyold + dzold*dzold); a += mass[type[i]]*(dxold*dxold + dyold*dyold + dzold*dzold);
b += mass[type[i]]*(dx *dxold + dy *dyold + dz *dzold); b += mass[type[i]]*(dx *dxold + dy *dyold + dz *dzold);
e += mass[type[i]]*(dx *dx + dy *dy + dz *dz); e += mass[type[i]]*(dx *dx + dy *dy + dz *dz);
@ -258,12 +244,10 @@ void FixTMD::initial_integrate(int vflag)
dxold = xold[i][0] - xf[i][0]; dxold = xold[i][0] - xf[i][0];
dyold = xold[i][1] - xf[i][1]; dyold = xold[i][1] - xf[i][1];
dzold = xold[i][2] - xf[i][2]; dzold = xold[i][2] - xf[i][2];
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; xback = unwrap[0] + gamma_back*dxold;
zbox = (image[i] >> IMG2BITS) - IMGMAX; yback = unwrap[1] + gamma_back*dyold;
xback = x[i][0] + xbox*xprd + gamma_back*dxold; zback = unwrap[2] + gamma_back*dzold;
yback = x[i][1] + ybox*yprd + gamma_back*dyold;
zback = x[i][2] + zbox*zprd + gamma_back*dzold;
dxkt = xback - xold[i][0]; dxkt = xback - xold[i][0];
dykt = yback - xold[i][1]; dykt = yback - xold[i][1];
dzkt = zback - xold[i][2]; dzkt = zback - xold[i][2];
@ -314,12 +298,7 @@ void FixTMD::initial_integrate(int vflag)
x[i][2] += gamma_forward*dzold; x[i][2] += gamma_forward*dzold;
v[i][2] += gamma_forward*dzold/dtv; v[i][2] += gamma_forward*dzold/dtv;
f[i][2] += gamma_forward*dzold/dtv/dtfm; f[i][2] += gamma_forward*dzold/dtv/dtfm;
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],xold[i]);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
xold[i][0] = x[i][0] + xbox*xprd;
xold[i][1] = x[i][1] + ybox*yprd;
xold[i][2] = x[i][2] + zbox*zprd;
} }
} }
} }

176
src/group.cpp
View File

@ -766,34 +766,27 @@ void Group::xcm(int igroup, double masstotal, double *cm)
double cmone[3]; double cmone[3];
cmone[0] = cmone[1] = cmone[2] = 0.0; cmone[0] = cmone[1] = cmone[2] = 0.0;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
double massone; double massone;
double unwrap[3];
if (rmass) { if (rmass) {
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = rmass[i]; massone = rmass[i];
cmone[0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
cmone[1] += (x[i][1] + ybox*yprd) * massone; cmone[0] += unwrap[0] * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone; cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
} }
} else { } else {
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = mass[type[i]]; massone = mass[type[i]];
cmone[0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
cmone[1] += (x[i][1] + ybox*yprd) * massone; cmone[0] += unwrap[0] * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone; cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
} }
} }
@ -827,34 +820,27 @@ void Group::xcm(int igroup, double masstotal, double *cm, int iregion)
double cmone[3]; double cmone[3];
cmone[0] = cmone[1] = cmone[2] = 0.0; cmone[0] = cmone[1] = cmone[2] = 0.0;
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
int xbox,ybox,zbox;
double massone; double massone;
double unwrap[3];
if (rmass) { if (rmass) {
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = rmass[i]; massone = rmass[i];
cmone[0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
cmone[1] += (x[i][1] + ybox*yprd) * massone; cmone[0] += unwrap[0] * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone; cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
} }
} else { } else {
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX;
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[i] >> IMG2BITS) - IMGMAX;
massone = mass[type[i]]; massone = mass[type[i]];
cmone[0] += (x[i][0] + xbox*xprd) * massone; domain->unmap(x[i],image[i],unwrap);
cmone[1] += (x[i][1] + ybox*yprd) * massone; cmone[0] += unwrap[0] * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone; cmone[1] += unwrap[1] * massone;
cmone[2] += unwrap[2] * massone;
} }
} }
@ -1102,21 +1088,16 @@ double Group::gyration(int igroup, double masstotal, double *cm)
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double rg = 0.0; double rg = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
rg += (dx*dx + dy*dy + dz*dz) * massone; rg += (dx*dx + dy*dy + dz*dz) * massone;
@ -1147,21 +1128,16 @@ double Group::gyration(int igroup, double masstotal, double *cm, int iregion)
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double rg = 0.0; double rg = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
rg += (dx*dx + dy*dy + dz*dz) * massone; rg += (dx*dx + dy*dy + dz*dz) * massone;
@ -1192,22 +1168,18 @@ void Group::angmom(int igroup, double *cm, double *lmom)
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double p[3]; double p[3];
p[0] = p[1] = p[2] = 0.0; p[0] = p[1] = p[2] = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
p[0] += massone * (dy*v[i][2] - dz*v[i][1]); p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
@ -1238,22 +1210,18 @@ void Group::angmom(int igroup, double *cm, double *lmom, int iregion)
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double p[3]; double p[3];
p[0] = p[1] = p[2] = 0.0; p[0] = p[1] = p[2] = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
p[0] += massone * (dy*v[i][2] - dz*v[i][1]); p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
@ -1280,22 +1248,18 @@ void Group::torque(int igroup, double *cm, double *tq)
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double tlocal[3]; double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[2] = 0.0; tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
tlocal[0] += dy*f[i][2] - dz*f[i][1]; tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2]; tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0]; tlocal[2] += dx*f[i][1] - dy*f[i][0];
@ -1321,22 +1285,18 @@ void Group::torque(int igroup, double *cm, double *tq, int iregion)
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double tlocal[3]; double tlocal[3];
tlocal[0] = tlocal[1] = tlocal[2] = 0.0; tlocal[0] = tlocal[1] = tlocal[2] = 0.0;
for (int i = 0; i < nlocal; i++) for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
tlocal[0] += dy*f[i][2] - dz*f[i][1]; tlocal[0] += dy*f[i][2] - dz*f[i][1];
tlocal[1] += dz*f[i][0] - dx*f[i][2]; tlocal[1] += dz*f[i][0] - dx*f[i][2];
tlocal[2] += dx*f[i][1] - dy*f[i][0]; tlocal[2] += dx*f[i][1] - dy*f[i][0];
@ -1364,11 +1324,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double ione[3][3]; double ione[3][3];
for (i = 0; i < 3; i++) for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++) for (j = 0; j < 3; j++)
@ -1376,12 +1334,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
ione[0][0] += massone * (dy*dy + dz*dz); ione[0][0] += massone * (dy*dy + dz*dz);
@ -1418,11 +1374,9 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
double *rmass = atom->rmass; double *rmass = atom->rmass;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz,massone; double dx,dy,dz,massone;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
double ione[3][3]; double ione[3][3];
for (i = 0; i < 3; i++) for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++) for (j = 0; j < 3; j++)
@ -1430,12 +1384,10 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion)
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0] - cm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - cm[1];
dx = (x[i][0] + xbox*xprd) - cm[0]; dz = unwrap[2] - cm[2];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (rmass) massone = rmass[i]; if (rmass) massone = rmass[i];
else massone = mass[type[i]]; else massone = mass[type[i]];
ione[0][0] += massone * (dy*dy + dz*dz); ione[0][0] += massone * (dy*dy + dz*dz);

15
src/velocity.cpp
View File

@ -718,20 +718,15 @@ void Velocity::zero_rotation()
tagint *image = atom->image; tagint *image = atom->image;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int xbox,ybox,zbox;
double dx,dy,dz; double dx,dy,dz;
double xprd = domain->xprd; double unwrap[3];
double yprd = domain->yprd;
double zprd = domain->zprd;
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
xbox = (image[i] & IMGMASK) - IMGMAX; domain->unmap(x[i],image[i],unwrap);
ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX; dx = unwrap[0]- xcm[0];
zbox = (image[i] >> IMG2BITS) - IMGMAX; dy = unwrap[1] - xcm[1];
dx = (x[i][0] + xbox*xprd) - xcm[0]; dz = unwrap[2] - xcm[2];
dy = (x[i][1] + ybox*yprd) - xcm[1];
dz = (x[i][2] + zbox*zprd) - xcm[2];
v[i][0] -= omega[1]*dz - omega[2]*dy; v[i][0] -= omega[1]*dz - omega[2]*dy;
v[i][1] -= omega[2]*dx - omega[0]*dz; v[i][1] -= omega[2]*dx - omega[0]*dz;
v[i][2] -= omega[0]*dy - omega[1]*dx; v[i][2] -= omega[0]*dy - omega[1]*dx;