jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2012-07-25 15:15:41 +00:00
parent 5e52200ad3
commit 912cfaca75
33 changed files with 649 additions and 651 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/atom.cpp
View File

@ -136,7 +136,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
map_style = 0;
map_tag_max = 0;
map_nhash = 0;
smax = 0;
sametag = NULL;
map_array = NULL;
@ -599,11 +599,11 @@ void Atom::data_atoms(int n, char *buf)
}
if (imageflag)
imagedata =
imagedata =
(((tagint) atoi(values[iptr+2]) + IMGMAX & IMGMASK) << IMG2BITS) |
(((tagint) atoi(values[iptr+1]) + IMGMAX & IMGMASK) << IMGBITS) |
(atoi(values[iptr]) + IMGMAX & IMGMASK);
else imagedata = ((tagint) IMGMAX << IMG2BITS) |
else imagedata = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
xdata[0] = atof(values[xptr]);

128
src/atom_vec_atomic.cpp
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
@ -46,7 +46,7 @@ AtomVecAtomic::AtomVecAtomic(LAMMPS *lmp, int narg, char **arg) :
/* ----------------------------------------------------------------------
grow atom arrays
n = 0 grows arrays by DELTA
n > 0 allocates arrays to size n
n > 0 allocates arrays to size n
------------------------------------------------------------------------- */
void AtomVecAtomic::grow(int n)
@ -66,7 +66,7 @@ void AtomVecAtomic::grow(int n)
f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
}
@ -99,14 +99,14 @@ void AtomVecAtomic::copy(int i, int j, int delflag)
v[j][2] = v[i][2];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j);
}
/* ---------------------------------------------------------------------- */
int AtomVecAtomic::pack_comm(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
@ -142,7 +142,7 @@ int AtomVecAtomic::pack_comm(int n, int *list, double *buf,
/* ---------------------------------------------------------------------- */
int AtomVecAtomic::pack_comm_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
@ -170,32 +170,32 @@ int AtomVecAtomic::pack_comm_vel(int n, int *list, double *buf,
}
if (!deform_vremap) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
}
}
}
@ -269,7 +269,7 @@ void AtomVecAtomic::unpack_reverse(int n, int *list, double *buf)
/* ---------------------------------------------------------------------- */
int AtomVecAtomic::pack_border(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz;
@ -311,7 +311,7 @@ int AtomVecAtomic::pack_border(int n, int *list, double *buf,
/* ---------------------------------------------------------------------- */
int AtomVecAtomic::pack_border_vel(int n, int *list, double *buf,
int pbc_flag, int *pbc)
int pbc_flag, int *pbc)
{
int i,j,m;
double dx,dy,dz,dvx,dvy,dvz;
@ -342,38 +342,38 @@ int AtomVecAtomic::pack_border_vel(int n, int *list, double *buf,
}
if (!deform_vremap) {
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
} else {
dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
dvz = pbc[2]*h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
j = list[i];
buf[m++] = x[j][0] + dx;
buf[m++] = x[j][1] + dy;
buf[m++] = x[j][2] + dz;
buf[m++] = tag[j];
buf[m++] = type[j];
buf[m++] = mask[j];
if (mask[i] & deform_groupbit) {
buf[m++] = v[j][0] + dvx;
buf[m++] = v[j][1] + dvy;
buf[m++] = v[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
}
}
}
}
@ -423,7 +423,7 @@ void AtomVecAtomic::unpack_border_vel(int n, int first, double *buf)
/* ----------------------------------------------------------------------
pack data for atom I for sending to another proc
xyz must be 1st 3 values, so comm::exchange() can test on them
xyz must be 1st 3 values, so comm::exchange() can test on them
------------------------------------------------------------------------- */
int AtomVecAtomic::pack_exchange(int i, double *buf)
@ -441,7 +441,7 @@ int AtomVecAtomic::pack_exchange(int i, double *buf)
*((tagint *) &buf[m++]) = image[i];
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
buf[0] = m;
@ -468,9 +468,9 @@ int AtomVecAtomic::unpack_exchange(double *buf)
image[nlocal] = *((tagint *) &buf[m++]);
if (atom->nextra_grow)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
m += modify->fix[atom->extra_grow[iextra]]->
unpack_exchange(nlocal,&buf[m]);
unpack_exchange(nlocal,&buf[m]);
atom->nlocal++;
return m;
@ -489,9 +489,9 @@ int AtomVecAtomic::size_restart()
int n = 11 * nlocal;
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
for (i = 0; i < nlocal; i++)
n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
return n;
}
@ -499,7 +499,7 @@ int AtomVecAtomic::size_restart()
/* ----------------------------------------------------------------------
pack atom I's data for restart file including extra quantities
xyz must be 1st 3 values, so that read_restart can test on them
molecular types may be negative, but write as positive
molecular types may be negative, but write as positive
------------------------------------------------------------------------- */
int AtomVecAtomic::pack_restart(int i, double *buf)
@ -517,7 +517,7 @@ int AtomVecAtomic::pack_restart(int i, double *buf)
buf[m++] = v[i][2];
if (atom->nextra_restart)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
buf[0] = m;
@ -575,7 +575,7 @@ void AtomVecAtomic::create_atom(int itype, double *coord)
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;
@ -617,7 +617,7 @@ void AtomVecAtomic::data_atom(double *coord, tagint imagetmp, char **values)
}
/* ----------------------------------------------------------------------
return # of bytes of allocated memory
return # of bytes of allocated memory
------------------------------------------------------------------------- */
bigint AtomVecAtomic::memory_usage()

2
src/atom_vec_charge.cpp
View File

@ -624,7 +624,7 @@ void AtomVecCharge::create_atom(int itype, double *coord)
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;

2
src/atom_vec_ellipsoid.cpp
View File

@ -1094,7 +1094,7 @@ void AtomVecEllipsoid::create_atom(int itype, double *coord)
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;

2
src/atom_vec_line.cpp
View File

@ -952,7 +952,7 @@ void AtomVecLine::create_atom(int itype, double *coord)
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;

2
src/atom_vec_sphere.cpp
View File

@ -922,7 +922,7 @@ void AtomVecSphere::create_atom(int itype, double *coord)
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;

2
src/atom_vec_tri.cpp
View File

@ -1303,7 +1303,7 @@ void AtomVecTri::create_atom(int itype, double *coord)
x[nlocal][1] = coord[1];
x[nlocal][2] = coord[2];
mask[nlocal] = 1;
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;

10
src/balance.cpp
View File

@ -707,7 +707,7 @@ int Balance::adjust(int n, double *split)
hi[i] = hi[i+1];
hisum[i] = hisum[i+1];
}
int change = 0;
for (int i = 1; i < n; i++)
if (sum[i] != target[i]) {
@ -888,11 +888,11 @@ void Balance::dumpout(bigint tstep, FILE *bfp)
fprintf(bfp,"%d\n",nprocs);
if (dimension == 2) fprintf(bfp,"ITEM: SQUARES\n");
else fprintf(bfp,"ITEM: CUBES\n");
int nx = comm->procgrid[0] + 1;
int ny = comm->procgrid[1] + 1;
int nz = comm->procgrid[2] + 1;
if (dimension == 2) {
int m = 0;
for (int j = 0; j < comm->procgrid[1]; j++)
@ -904,7 +904,7 @@ void Balance::dumpout(bigint tstep, FILE *bfp)
fprintf(bfp,"%d %d %d %d %d %d\n",m+1,m+1,c1,c2,c3,c4);
m++;
}
} else {
int m = 0;
for (int k = 0; k < comm->procgrid[2]; k++)
@ -1000,7 +1000,7 @@ void Balance::debug_output(int idim, int m, int np, double *split)
for (i = 0; i <= np; i++) printf(" " BIGINT_FORMAT,target[i]);
printf("\n");
printf(" Actual cut:");
for (i = 0; i <= np; i++)
for (i = 0; i <= np; i++)
printf(" %g",boxlo[bdim[idim]] + split[i]*prd[bdim[idim]]);
printf("\n");
printf(" Split:");

2
src/balance.h
View File

@ -47,7 +47,7 @@ class Balance : protected Pointers {
char bstr[4];
int ndim; // length of balance string bstr
int *bdim; // XYZ for each character in bstr
int *bdim; // XYZ for each character in bstr
bigint *count; // counts for slices in one dim
bigint *onecount; // work vector of counts in one dim
bigint *sum; // cummulative count for slices in one dim

2
src/compute_gyration.cpp
View File

@ -88,7 +88,7 @@ void ComputeGyration::compute_vector()
double xprd = domain->xprd;
double yprd = domain->yprd;
double zprd = domain->zprd;
double rg[6];
rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0;

6
src/compute_gyration_molecule.cpp
View File

@ -179,7 +179,7 @@ void ComputeGyrationMolecule::compute_array()
molcom();
for (i = 0; i < nmolecules; i++)
for (j = 0; j < 6; j++)
for (j = 0; j < 6; j++)
rgt[i][j] = 0.0;
double **x = atom->x;
@ -220,8 +220,8 @@ void ComputeGyrationMolecule::compute_array()
MPI_Allreduce(&rgt[0][0],&array[0][0],nmolecules*6,
MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < nmolecules; i++)
for (j = 0; j < 6; j++)
for (i = 0; i < nmolecules; i++)
for (j = 0; j < 6; j++)
array[i][j] /= masstotal[i];
}

2
src/create_atoms.cpp
View File

@ -271,7 +271,7 @@ void CreateAtoms::add_single()
// remap atom if requested
if (remapflag) {
tagint imagetmp = ((tagint) IMGMAX << IMG2BITS) |
tagint imagetmp = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(xone,imagetmp);
}

6
src/domain.cpp
View File

@ -535,7 +535,7 @@ void Domain::box_too_small_check()
// NOTE: if box is tiny (less than 2 * bond-length),
// the check itself may compute bad bond lengths
// not sure how to account for that extreme case
int *num_bond = atom->num_bond;
int **bond_atom = atom->bond_atom;
double **x = atom->x;
@ -724,7 +724,7 @@ int Domain::closest_image(int i, int j)
double delz = xi[2] - x[j][2];
double rsqmin = delx*delx + dely*dely + delz*delz;
double rsq;
while (sametag[j] >= 0) {
j = sametag[j];
delx = xi[0] - x[j][0];
@ -1138,7 +1138,7 @@ void Domain::image_flip(int m, int n, int p)
xbox -= m*ybox + n*zbox;
image[i] = ((zbox + (tagint) IMGMAX & IMGMASK) << IMG2BITS) |
((ybox + (tagint) IMGMAX & IMGMASK) << IMGBITS) |
((ybox + (tagint) IMGMAX & IMGMASK) << IMGBITS) |
(xbox + IMGMAX & IMGMASK);
}
}

8
src/fix_aveforce.cpp
View File

@ -128,21 +128,21 @@ void FixAveForce::init()
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0)
if (xvar < 0)
error->all(FLERR,"Variable name for fix aveforce does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else error->all(FLERR,"Variable for fix aveforce is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0)
if (yvar < 0)
error->all(FLERR,"Variable name for fix aveforce does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else error->all(FLERR,"Variable for fix aveforce is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0)
if (zvar < 0)
error->all(FLERR,"Variable name for fix aveforce does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else error->all(FLERR,"Variable for fix aveforce is invalid style");
@ -152,7 +152,7 @@ void FixAveForce::init()
if (iregion >= 0) {
iregion = domain->find_region(idregion);
if (iregion == -1)
if (iregion == -1)
error->all(FLERR,"Region ID for fix aveforce does not exist");
}

12
src/fix_balance.cpp
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
@ -52,13 +52,13 @@ FixBalance::FixBalance(LAMMPS *lmp, int narg, char **arg) :
nitermax = atoi(arg[5]);
thresh = atof(arg[6]);
if (nevery < 0 || nitermax <= 0 || thresh < 1.0)
if (nevery < 0 || nitermax <= 0 || thresh < 1.0)
error->all(FLERR,"Illegal fix balance command");
for (int i = 0; i < strlen(bstr); i++) {
if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z')
if (bstr[i] != 'x' && bstr[i] != 'y' && bstr[i] != 'z')
error->all(FLERR,"Fix balance string is invalid");
if (bstr[i] == 'z' && dimension == 2)
if (bstr[i] == 'z' && dimension == 2)
error->all(FLERR,"Fix balance string is invalid for 2d simulation");
for (int j = i+1; j < strlen(bstr); j++)
if (bstr[i] == bstr[j])
@ -220,7 +220,7 @@ void FixBalance::rebalance()
itercount = balance->dynamic();
// output of final result
if (fp) balance->dumpout(update->ntimestep,fp);
// reset comm->uniform flag
@ -236,7 +236,7 @@ void FixBalance::rebalance()
// move atoms to new processors via irregular()
// only needed if migrate_check() says an atom moves to far,
// else allow caller's comm->exchange() to do it
if (domain->triclinic) domain->x2lamda(atom->nlocal);
if (irregular->migrate_check()) irregular->migrate_atoms();
if (domain->triclinic) domain->lamda2x(atom->nlocal);

2
src/fix_balance.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.

2
src/fix_nh.cpp
View File

@ -438,7 +438,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
// reneighboring only forced if flips will occur due to shape changes
if (p_flag[3] || p_flag[4] || p_flag[5]) force_reneighbor = 1;
if (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0)
if (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0)
force_reneighbor = 1;
// convert input periods to frequencies

16
src/fix_rigid.cpp
View File

@ -1651,7 +1651,7 @@ void FixRigid::setup_bodies()
// then remap the xcm of each body back into simulation box if needed
for (ibody = 0; ibody < nbody; ibody++)
imagebody[ibody] = ((tagint) IMGMAX << IMG2BITS) |
imagebody[ibody] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
pre_neighbor();
@ -2003,7 +2003,7 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
FILE *fp;
char *eof,*start,*next,*buf;
char line[MAXLINE];
if (me == 0) {
fp = fopen(infile,"r");
if (fp == NULL) {
@ -2052,7 +2052,7 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
if (nwords != ATTRIBUTE_PERBODY)
error->all(FLERR,"Incorrect rigid body format in fix rigid file");
// loop over lines of rigid body attributes
// tokenize the line into values
// id = rigid body ID
@ -2062,19 +2062,19 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
for (int i = 0; i < nchunk; i++) {
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
for (j = 1; j < nwords; j++)
values[j] = strtok(NULL," \t\n\r\f");
id = atoi(values[0]);
if (rstyle == MOLECULE) {
if (id <= 0 || id > maxmol)
if (id <= 0 || id > maxmol)
error->all(FLERR,"Invalid rigid body ID in fix rigid file");
id = mol2body[id];
} else id--;
if (id < 0 || id >= nbody)
if (id < 0 || id >= nbody)
error->all(FLERR,"Invalid rigid body ID in fix rigid file");
inbody[id] = 1;
@ -2094,7 +2094,7 @@ void FixRigid::readfile(int which, double *vec, double **array, int *inbody)
buf = next + 1;
}
nread += nchunk;
}

2
src/fix_shake.cpp
View File

@ -126,7 +126,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
} else if (mode == 'm') {
double massone = atof(arg[next]);
if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake");
if (nmass == atom->ntypes)
if (nmass == atom->ntypes)
error->all(FLERR,"Too many masses for fix shake");
mass_list[nmass++] = massone;

2
src/lmptype.h
View File

@ -5,7 +5,7 @@
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.

2
src/neigh_half_bin.cpp
View File

@ -203,7 +203,7 @@ void Neighbor::half_bin_no_newton_ghost(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
which = find_special(special[i],nspecial[i],tag[j]);

10
src/neigh_half_nsq.cpp
View File

@ -179,12 +179,12 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
if (includegroup && !(mask[j] & bitmask)) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) {
if (molecular) {
which = find_special(special[i],nspecial[i],tag[j]);
@ -201,16 +201,16 @@ void Neighbor::half_nsq_no_newton_ghost(NeighList *list)
if (includegroup && !(mask[j] & bitmask)) continue;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= cutneighsq[itype][jtype]) neighptr[n++] = j;
}
}
ilist[inum++] = i;
firstneigh[i] = neighptr;
numneigh[i] = n;

6
src/neighbor.cpp
View File

@ -493,7 +493,7 @@ void Neighbor::init()
lists[i]->listskip = lists[requests[i]->otherlist];
lists[i]->copy_skip_info(requests[i]->iskip,requests[i]->ijskip);
} else if (requests[i]->half_from_full)
} else if (requests[i]->half_from_full)
lists[i]->listfull = lists[i-1];
else if (requests[i]->granhistory) {
@ -852,7 +852,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
else pb = &Neighbor::half_bin_no_newton_ghost;
} else if (triclinic == 0) {
pb = &Neighbor::half_bin_newton;
} else if (triclinic == 1)
} else if (triclinic == 1)
pb = &Neighbor::half_bin_newton_tri;
} else if (rq->newton == 1) {
if (triclinic == 0) pb = &Neighbor::half_bin_newton;
@ -966,7 +966,7 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
else pb = &Neighbor::half_bin_no_newton_ghost_omp;
} else if (triclinic == 0) {
pb = &Neighbor::half_bin_newton_omp;
} else if (triclinic == 1)
} else if (triclinic == 1)
pb = &Neighbor::half_bin_newton_tri_omp;
} else if (rq->newton == 1) {
if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp;

8
src/output.cpp
View File

@ -472,7 +472,7 @@ void Output::reset_timestep(bigint ntimestep)
if (restart_every_single) {
next_restart_single =
(ntimestep/restart_every_single)*restart_every_single;
if (next_restart_single < ntimestep)
if (next_restart_single < ntimestep)
next_restart_single += restart_every_single;
} else {
modify->clearstep_compute();
@ -491,7 +491,7 @@ void Output::reset_timestep(bigint ntimestep)
if (restart_every_double) {
next_restart_double =
(ntimestep/restart_every_double)*restart_every_double;
if (next_restart_double < ntimestep)
if (next_restart_double < ntimestep)
next_restart_double += restart_every_double;
} else {
modify->clearstep_compute();
@ -507,7 +507,7 @@ void Output::reset_timestep(bigint ntimestep)
} else next_restart_double = update->laststep + 1;
next_restart = MIN(next_restart_single,next_restart_double);
if (var_thermo) {
modify->clearstep_compute();
update->ntimestep--;
@ -522,7 +522,7 @@ void Output::reset_timestep(bigint ntimestep)
next_thermo = (ntimestep/thermo_every)*thermo_every;
if (next_thermo < ntimestep) next_thermo += thermo_every;
next_thermo = MIN(next_thermo,update->laststep);
} else next_thermo = update->laststep;
} else next_thermo = update->laststep;
next = MIN(next_dump_any,next_restart);
next = MIN(next,next_thermo);

2
src/pair.h
View File

@ -57,7 +57,7 @@ class Pair : protected Pointers {
int evflag; // energy,virial settings
int eflag_either,eflag_global,eflag_atom;
int vflag_either,vflag_global,vflag_atom;
int ncoultablebits; // size of Coulomb table, accessed by KSpace
int nextra; // # of extra quantities pair style calculates

885
src/pair_lj_cut_coul_cut.cpp
View File

@ -1,443 +1,442 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_cut_coul_cut.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairLJCutCoulCut::PairLJCutCoulCut(LAMMPS *lmp) : Pair(lmp) {}
/* ---------------------------------------------------------------------- */
PairLJCutCoulCut::~PairLJCutCoulCut()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut_lj);
memory->destroy(cut_ljsq);
memory->destroy(cut_coul);
memory->destroy(cut_coulsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
}
}
/* ---------------------------------------------------------------------- */
void PairLJCutCoulCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
if (rsq < cut_coulsq[itype][jtype])
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
if (rsq < cut_coulsq[itype][jtype])
ecoul = factor_coul * qqrd2e * qtmp*q[j]*sqrt(r2inv);
else ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJCutCoulCut::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
memory->create(cut_coul,n+1,n+1,"pair:cut_coul");
memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJCutCoulCut::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
cut_lj_global = force->numeric(arg[0]);
if (narg == 1) cut_coul_global = cut_lj_global;
else cut_coul_global = force->numeric(arg[1]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) {
cut_lj[i][j] = cut_lj_global;
cut_coul[i][j] = cut_coul_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJCutCoulCut::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double epsilon_one = force->numeric(arg[2]);
double sigma_one = force->numeric(arg[3]);
double cut_lj_one = cut_lj_global;
double cut_coul_one = cut_coul_global;
if (narg >= 5) cut_coul_one = cut_lj_one = force->numeric(arg[4]);
if (narg == 6) cut_coul_one = force->numeric(arg[5]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut_lj[i][j] = cut_lj_one;
cut_coul[i][j] = cut_coul_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJCutCoulCut::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/cut/coul/cut requires atom attribute q");
neighbor->request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJCutCoulCut::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]);
}
double cut = MAX(cut_lj[i][j],cut_coul[i][j]);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j];
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
} else offset[i][j] = 0.0;
cut_ljsq[j][i] = cut_ljsq[i][j];
cut_coulsq[j][i] = cut_coulsq[i][j];
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
offset[j][i] = offset[i][j];
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
if (tail_flag) {
int *type = atom->type;
int nlocal = atom->nlocal;
double count[2],all[2];
count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0;
}
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3;
double rc9 = rc3*rc6;
etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
}
return cut;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCutCoulCut::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
fwrite(&cut_coul[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCutCoulCut::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut_lj[i][j],sizeof(double),1,fp);
fread(&cut_coul[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCutCoulCut::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCutCoulCut::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairLJCutCoulCut::single(int i, int j, int itype, int jtype,
double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r6inv,forcecoul,forcelj,phicoul,philj;
r2inv = 1.0/rsq;
if (rsq < cut_coulsq[itype][jtype])
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
double eng = 0.0;
if (rsq < cut_coulsq[itype][jtype]) {
phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
eng += factor_coul*phicoul;
}
if (rsq < cut_ljsq[itype][jtype]) {
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
eng += factor_lj*philj;
}
return eng;
}
/* ---------------------------------------------------------------------- */
void *PairLJCutCoulCut::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
return NULL;
}
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_lj_cut_coul_cut.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairLJCutCoulCut::PairLJCutCoulCut(LAMMPS *lmp) : Pair(lmp) {}
/* ---------------------------------------------------------------------- */
PairLJCutCoulCut::~PairLJCutCoulCut()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut_lj);
memory->destroy(cut_ljsq);
memory->destroy(cut_coul);
memory->destroy(cut_coulsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
memory->destroy(offset);
}
}
/* ---------------------------------------------------------------------- */
void PairLJCutCoulCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double qqrd2e = force->qqrd2e;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
qtmp = q[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
if (rsq < cut_coulsq[itype][jtype])
forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
if (rsq < cut_coulsq[itype][jtype])
ecoul = factor_coul * qqrd2e * qtmp*q[j]*sqrt(r2inv);
else ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,ecoul,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairLJCutCoulCut::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
memory->create(cut_coul,n+1,n+1,"pair:cut_coul");
memory->create(cut_coulsq,n+1,n+1,"pair:cut_coulsq");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairLJCutCoulCut::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
cut_lj_global = force->numeric(arg[0]);
if (narg == 1) cut_coul_global = cut_lj_global;
else cut_coul_global = force->numeric(arg[1]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i+1; j <= atom->ntypes; j++)
if (setflag[i][j]) {
cut_lj[i][j] = cut_lj_global;
cut_coul[i][j] = cut_coul_global;
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairLJCutCoulCut::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(arg[0],atom->ntypes,ilo,ihi);
force->bounds(arg[1],atom->ntypes,jlo,jhi);
double epsilon_one = force->numeric(arg[2]);
double sigma_one = force->numeric(arg[3]);
double cut_lj_one = cut_lj_global;
double cut_coul_one = cut_coul_global;
if (narg >= 5) cut_coul_one = cut_lj_one = force->numeric(arg[4]);
if (narg == 6) cut_coul_one = force->numeric(arg[5]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut_lj[i][j] = cut_lj_one;
cut_coul[i][j] = cut_coul_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLJCutCoulCut::init_style()
{
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/cut/coul/cut requires atom attribute q");
neighbor->request(this);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairLJCutCoulCut::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
cut_coul[i][j] = mix_distance(cut_coul[i][i],cut_coul[j][j]);
}
double cut = MAX(cut_lj[i][j],cut_coul[i][j]);
cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];
cut_coulsq[i][j] = cut_coul[i][j] * cut_coul[i][j];
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
if (offset_flag) {
double ratio = sigma[i][j] / cut_lj[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
} else offset[i][j] = 0.0;
cut_ljsq[j][i] = cut_ljsq[i][j];
cut_coulsq[j][i] = cut_coulsq[i][j];
lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j];
lj3[j][i] = lj3[i][j];
lj4[j][i] = lj4[i][j];
offset[j][i] = offset[i][j];
// compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce
if (tail_flag) {
int *type = atom->type;
int nlocal = atom->nlocal;
double count[2],all[2];
count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0;
}
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
double sig2 = sigma[i][j]*sigma[i][j];
double sig6 = sig2*sig2*sig2;
double rc3 = cut_lj[i][j]*cut_lj[i][j]*cut_lj[i][j];
double rc6 = rc3*rc3;
double rc9 = rc3*rc6;
etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9);
ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] *
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
}
return cut;
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCutCoulCut::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp);
fwrite(&sigma[i][j],sizeof(double),1,fp);
fwrite(&cut_lj[i][j],sizeof(double),1,fp);
fwrite(&cut_coul[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCutCoulCut::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&epsilon[i][j],sizeof(double),1,fp);
fread(&sigma[i][j],sizeof(double),1,fp);
fread(&cut_lj[i][j],sizeof(double),1,fp);
fread(&cut_coul[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairLJCutCoulCut::write_restart_settings(FILE *fp)
{
fwrite(&cut_lj_global,sizeof(double),1,fp);
fwrite(&cut_coul_global,sizeof(double),1,fp);
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairLJCutCoulCut::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_lj_global,sizeof(double),1,fp);
fread(&cut_coul_global,sizeof(double),1,fp);
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
double PairLJCutCoulCut::single(int i, int j, int itype, int jtype,
double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,r6inv,forcecoul,forcelj,phicoul,philj;
r2inv = 1.0/rsq;
if (rsq < cut_coulsq[itype][jtype])
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
else forcecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
} else forcelj = 0.0;
fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;
double eng = 0.0;
if (rsq < cut_coulsq[itype][jtype]) {
phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
eng += factor_coul*phicoul;
}
if (rsq < cut_ljsq[itype][jtype]) {
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
eng += factor_lj*philj;
}
return eng;
}
/* ---------------------------------------------------------------------- */
void *PairLJCutCoulCut::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
return NULL;
}

148
src/pair_lj_cut_coul_cut.h
View File

@ -1,74 +1,74 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/cut/coul/cut,PairLJCutCoulCut)
#else
#ifndef LMP_PAIR_LJ_CUT_COUL_CUT_H
#define LMP_PAIR_LJ_CUT_COUL_CUT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCutCoulCut : public Pair {
public:
PairLJCutCoulCut(class LAMMPS *);
virtual ~PairLJCutCoulCut();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
void *extract(const char *, int &);
protected:
double cut_lj_global,cut_coul_global;
double **cut_lj,**cut_ljsq;
double **cut_coul,**cut_coulsq;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style lj/cut/coul/cut requires atom attribute q
The atom style defined does not have this attribute.
*/
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(lj/cut/coul/cut,PairLJCutCoulCut)
#else
#ifndef LMP_PAIR_LJ_CUT_COUL_CUT_H
#define LMP_PAIR_LJ_CUT_COUL_CUT_H
#include "pair.h"
namespace LAMMPS_NS {
class PairLJCutCoulCut : public Pair {
public:
PairLJCutCoulCut(class LAMMPS *);
virtual ~PairLJCutCoulCut();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
virtual double single(int, int, int, int, double, double, double, double &);
void *extract(const char *, int &);
protected:
double cut_lj_global,cut_coul_global;
double **cut_lj,**cut_ljsq;
double **cut_coul,**cut_coulsq;
double **epsilon,**sigma;
double **lj1,**lj2,**lj3,**lj4,**offset;
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style lj/cut/coul/cut requires atom attribute q
The atom style defined does not have this attribute.
*/

5
src/read_dump.cpp
View File

@ -729,7 +729,7 @@ void ReadDump::process_atoms(int n)
// replace image flag in case changed by ix,iy,iz fields
image[m] = ((tagint) xbox << IMG2BITS) |
image[m] = ((tagint) xbox << IMG2BITS) |
((tagint) ybox << IMGBITS) | zbox;
}
}
@ -814,7 +814,7 @@ void ReadDump::process_atoms(int n)
// replace image flag in case changed by ix,iy,iz fields
image[m] = ((tagint) xbox << IMG2BITS) |
image[m] = ((tagint) xbox << IMG2BITS) |
((tagint) ybox << IMGBITS) | zbox;
}
}
@ -881,4 +881,3 @@ double ReadDump::zfield(int i, int j)
if (scaled == UNSCALED) return fields[i][j];
return fields[i][j]*zprd + zlo;
}

4
src/reader_xyz.cpp
View File

@ -127,7 +127,7 @@ bigint ReaderXYZ::read_header(double box[3][3], int &triclinic,
memory->create(fieldindex,nfield,"read_dump:fieldindex");
// for xyz we know nothing about the style of coordinates,
// for xyz we know nothing about the style of coordinates,
// so the caller has to set the proper flag.
xflag = scaledflag;
@ -148,7 +148,7 @@ bigint ReaderXYZ::read_header(double box[3][3], int &triclinic,
fieldflag = 1;
}
}
return natoms;
}

2
src/replicate.cpp
View File

@ -278,7 +278,7 @@ void Replicate::command(int narg, char **arg)
m = 0;
while (m < n) {
image = ((tagint) IMGMAX << IMG2BITS) |
image = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
if (triclinic == 0) {
x[0] = buf[m+1] + ix*old_xprd;

4
src/rerun.cpp
View File

@ -108,7 +108,7 @@ void Rerun::command(int narg, char **arg)
int nremain = narg - iarg - 1;
if (nremain <= 0) error->all(FLERR,"Illegal rerun command");
if (first > last) error->all(FLERR,"Illegal rerun command");
if (startflag && stopflag && start > stop)
if (startflag && stopflag && start > stop)
error->all(FLERR,"Illegal rerun command");
// pass list of filenames to ReadDump
@ -117,7 +117,7 @@ void Rerun::command(int narg, char **arg)
ReadDump *rd = new ReadDump(lmp);
rd->store_files(nfile,arg);
if (nremain)
if (nremain)
nremain = rd->fields_and_keywords(nremain,&arg[narg-nremain]);
else nremain = rd->fields_and_keywords(0,NULL);
if (nremain) rd->setup_reader(nremain,&arg[narg-nremain]);

2
src/set.cpp
View File

@ -497,7 +497,7 @@ void Set::set(int keyword)
if (yimageflag) ybox = yimage;
if (zimageflag) zbox = zimage;
atom->image[i] = ((zbox + (tagint) IMGMAX & IMGMASK) << IMG2BITS) |
((ybox + (tagint) IMGMAX & IMGMASK) << IMGBITS) |
((ybox + (tagint) IMGMAX & IMGMASK) << IMGBITS) |
(xbox + IMGMAX & IMGMASK);
// set dipole moment

6
src/variable.cpp
View File

@ -164,7 +164,7 @@ void Variable::set(int narg, char **arg)
} else if (narg == 4 || (narg == 5 && strcmp(arg[4],"pad") == 0)) {
nfirst = atoi(arg[2]);
nlast = atoi(arg[3]);
if (nfirst > nlast || nlast < 0)
if (nfirst > nlast || nlast < 0)
error->all(FLERR,"Illegal variable command");
if (narg == 5 && strcmp(arg[4],"pad") == 0) {
char digits[12];
@ -924,7 +924,7 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->type = ATOMARRAY;
if (compute->array_atom)
newtree->array = &compute->array_atom[0][index1-1];
else
else
newtree->array = NULL;
newtree->nstride = compute->size_peratom_cols;
newtree->left = newtree->middle = newtree->right = NULL;
@ -1093,7 +1093,7 @@ double Variable::evaluate(char *str, Tree **tree)
newtree->type = ATOMARRAY;
if (fix->array_atom)
newtree->array = &fix->array_atom[0][index1-1];
else
else
newtree->array = NULL;
newtree->nstride = fix->size_peratom_cols;
newtree->left = newtree->middle = newtree->right = NULL;