lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2009-02-12 21:36:52 +00:00
parent 0b60bfbf51
commit 94ddc19477
39 changed files with 1506 additions and 1426 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/COLLOID/pair_colloid.cpp
View File

@ -521,11 +521,3 @@ double PairColloid::single(int i, int j, int itype, int jtype, double rsq,
return factor_lj*phi;
}
/* ---------------------------------------------------------------------- */
void *PairColloid::extract(char *str)
{
if (strcmp(str,"diameter") == 0) return (void *) diameter;
return NULL;
}

1
src/COLLOID/pair_colloid.h
View File

@ -31,7 +31,6 @@ class PairColloid : public Pair {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double single(int, int, int, int, double, double, double, double &);
void *extract(char *);
private:
double cut_global;

2
src/GRANULAR/atom_vec_granular.cpp
View File

@ -604,6 +604,7 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values)
density[nlocal] = atof(values[3]);
rmass[nlocal] = 4.0*PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
if (rmass[nlocal] <= 0.0) error->one("Invalid mass value");
x[nlocal][0] = coord[0];
x[nlocal][1] = coord[1];
@ -633,6 +634,7 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values)
density[nlocal] = atof(values[1]);
rmass[nlocal] = 4.0*PI/3.0 *
radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];
if (rmass[nlocal] <= 0.0) error->one("Invalid mass value");
v[nlocal][0] = 0.0;
v[nlocal][1] = 0.0;

12
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
View File

@ -54,6 +54,10 @@ PairLJCutCoulLongTIP4P::PairLJCutCoulLongTIP4P(LAMMPS *lmp) :
PairLJCutCoulLong(lmp)
{
single_enable = 0;
// TIP4P cannot compute virial as F dot r
// due to find_M() finding bonded H atoms which are not near O atom
no_virial_compute = 1;
}
@ -76,15 +80,9 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
float rsq;
int *int_rsq = (int *) &rsq;
// if global component of incoming vflag = 2, reset vflag as if it were 1
// necessary since TIP4P cannot compute virial as F dot r
// due to find_M() finding bonded H atoms which are not near O atom
if (vflag % 4 == 2) vflag = 1 + vflag/4 * 4;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = 0;
else evflag = vflag_fdotr = 0;
double **f = atom->f;
double **x = atom->x;

5
src/MANYBODY/pair_sw.cpp
View File

@ -122,9 +122,8 @@ void PairSW::compute(int eflag, int vflag)
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
else if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)
continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
jtype = map[type[j]];

5
src/MANYBODY/pair_tersoff.cpp
View File

@ -125,9 +125,8 @@ void PairTersoff::compute(int eflag, int vflag)
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < x[i][2]) continue;
else if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)
continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
jtype = map[type[j]];

1220
src/PERI/atom_vec_peri.cpp
View File

File diff suppressed because it is too large Load Diff

1042
src/PERI/pair_peri_pmb.cpp
View File

File diff suppressed because it is too large Load Diff

1
src/atom.h
View File

@ -73,6 +73,7 @@ class Atom : protected Pointers {
int **improper_atom1,**improper_atom2,**improper_atom3,**improper_atom4;
// per-atom array existence flags
// these can be checked before array is allocated
// customize by adding new flag
int molecule_flag;

2
src/compute_centro_atom.h
View File

@ -31,8 +31,8 @@ class ComputeCentroAtom : public Compute {
int nmax,maxneigh;
double *distsq;
int *nearest;
double *centro;
class NeighList *list;
double *centro;
void select(int, int, double *);
void select2(int, int, double *, int *);

9
src/compute_coord_atom.cpp
View File

@ -11,6 +11,7 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "compute_coord_atom.h"
@ -35,7 +36,8 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 4) error->all("Illegal compute coord/atom command");
cutoff = atof(arg[3]);
double cutoff = atof(arg[3]);
cutsq = cutoff*cutoff;
peratom_flag = 1;
size_peratom = 0;
@ -55,7 +57,9 @@ ComputeCoordAtom::~ComputeCoordAtom()
void ComputeCoordAtom::init()
{
if (force->pair == NULL || cutoff > force->pair->cutforce)
if (force->pair == NULL)
error->all("Compute coord/atom requires a pair style be defined");
if (sqrt(cutsq) > force->pair->cutforce)
error->all("Compute coord/atom cutoff is longer than pairwise cutoff");
// need an occasional full neighbor list
@ -116,7 +120,6 @@ void ComputeCoordAtom::compute_peratom()
double **x = atom->x;
int *mask = atom->mask;
int nall = atom->nlocal + atom->nghost;
double cutsq = cutoff*cutoff;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];

4
src/compute_coord_atom.h
View File

@ -29,9 +29,9 @@ class ComputeCoordAtom : public Compute {
private:
int nmax;
double cutoff;
double *coordination;
double cutsq;
class NeighList *list;
double *coordination;
};
}

16
src/compute_erotate_sphere.cpp
View File

@ -55,8 +55,8 @@ void ComputeERotateSphere::init()
if (atom->mass && !atom->shape)
error->all("Compute erotate/sphere requires atom attribute shape");
if (!atom->mass && (!atom->radius_flag || !atom->rmass_flag))
error->all("Compute erotate/sphere requires atom attributes radius, rmass");
if (atom->rmass && !atom->radius_flag)
error->all("Compute erotate/sphere requires atom attribute radius");
if (atom->mass) {
double *mass = atom->mass;
@ -86,16 +86,16 @@ double ComputeERotateSphere::compute_scalar()
double erotate = 0.0;
if (mass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * inertia[type[i]];
} else {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * inertia[type[i]];
}
MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);

11
src/compute_ke.cpp
View File

@ -55,16 +55,15 @@ double ComputeKE::compute_scalar()
double ke = 0.0;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
ke += mass[type[i]] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
ke += rmass[i] * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
}
MPI_Allreduce(&ke,&scalar,1,MPI_DOUBLE,MPI_SUM,world);

6
src/compute_ke_atom.cpp
View File

@ -80,10 +80,10 @@ void ComputeKEAtom::compute_peratom()
int *type = atom->type;
int nlocal = atom->nlocal;
if (mass)
if (rmass)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ke[i] = 0.5 * mvv2e * mass[type[i]] *
ke[i] = 0.5 * mvv2e * rmass[i] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
} else ke[i] = 0.0;
}
@ -91,7 +91,7 @@ void ComputeKEAtom::compute_peratom()
else
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
ke[i] = 0.5 * mvv2e * rmass[i] *
ke[i] = 0.5 * mvv2e * mass[type[i]] *
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
} else ke[i] = 0.0;
}

6
src/compute_stress_atom.cpp
View File

@ -196,10 +196,10 @@ void ComputeStressAtom::compute_peratom()
int *type = atom->type;
double mvv2e = force->mvv2e;
if (mass) {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
onemass = mvv2e * mass[type[i]];
onemass = mvv2e * rmass[i];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];
@ -211,7 +211,7 @@ void ComputeStressAtom::compute_peratom()
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
onemass = mvv2e * rmass[i];
onemass = mvv2e * mass[type[i]];
stress[i][0] += onemass*v[i][0]*v[i][0];
stress[i][1] += onemass*v[i][1]*v[i][1];
stress[i][2] += onemass*v[i][2]*v[i][2];

14
src/compute_temp.cpp
View File

@ -84,15 +84,15 @@ double ComputeTemp::compute_scalar()
double t = 0.0;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
}
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
@ -121,8 +121,8 @@ void ComputeTemp::compute_vector()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * v[i][0]*v[i][0];
t[1] += massone * v[i][1]*v[i][1];
t[2] += massone * v[i][2]*v[i][2];

12
src/compute_temp_com.cpp
View File

@ -102,12 +102,12 @@ double ComputeTempCOM::compute_scalar()
vthermal[0] = v[i][0] - vbias[0];
vthermal[1] = v[i][1] - vbias[1];
vthermal[2] = v[i][2] - vbias[2];
if (mass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
else
if (rmass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * rmass[i];
else
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
}
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
@ -144,8 +144,8 @@ void ComputeTempCOM::compute_vector()
vthermal[1] = v[i][1] - vbias[1];
vthermal[2] = v[i][2] - vbias[2];
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0];
t[1] += massone * vthermal[1]*vthermal[1];
t[2] += massone * vthermal[2]*vthermal[2];

15
src/compute_temp_deform.cpp
View File

@ -78,7 +78,8 @@ void ComputeTempDeform::init()
if (strcmp(modify->fix[i]->style,"deform") == 0) {
if (((FixDeform *) modify->fix[i])->remapflag == X_REMAP &&
comm->me == 0)
error->warning("Using compute temp/deform with inconsistent fix deform remap option");
error->warning("Using compute temp/deform with inconsistent "
"fix deform remap option");
break;
}
if (i == modify->nfix && comm->me == 0)
@ -131,12 +132,12 @@ double ComputeTempDeform::compute_scalar()
vthermal[0] = v[i][0] - vstream[0];
vthermal[1] = v[i][1] - vstream[1];
vthermal[2] = v[i][2] - vstream[2];
if (mass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
else
if (rmass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * rmass[i];
else
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
}
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
@ -178,8 +179,8 @@ void ComputeTempDeform::compute_vector()
vthermal[1] = v[i][1] - vstream[1];
vthermal[2] = v[i][2] - vstream[2];
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0];
t[1] += massone * vthermal[1]*vthermal[1];
t[2] += massone * vthermal[2]*vthermal[2];

16
src/compute_temp_partial.cpp
View File

@ -104,16 +104,16 @@ double ComputeTempPartial::compute_scalar()
double t = 0.0;
if (mass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] +
zflag*v[i][2]*v[i][2]) * mass[type[i]];
} else {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] +
zflag*v[i][2]*v[i][2]) * rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
t += (xflag*v[i][0]*v[i][0] + yflag*v[i][1]*v[i][1] +
zflag*v[i][2]*v[i][2]) * mass[type[i]];
}
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
@ -142,8 +142,8 @@ void ComputeTempPartial::compute_vector()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * xflag*v[i][0]*v[i][0];
t[1] += massone * yflag*v[i][1]*v[i][1];
t[2] += massone * zflag*v[i][2]*v[i][2];

12
src/compute_temp_ramp.cpp
View File

@ -170,12 +170,12 @@ double ComputeTempRamp::compute_scalar()
vthermal[1] = v[i][1];
vthermal[2] = v[i][2];
vthermal[v_dim] -= vramp;
if (mass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
else
if (rmass)
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * rmass[i];
else
t += (vthermal[0]*vthermal[0] + vthermal[1]*vthermal[1] +
vthermal[2]*vthermal[2]) * mass[type[i]];
}
MPI_Allreduce(&t,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
@ -215,8 +215,8 @@ void ComputeTempRamp::compute_vector()
vthermal[2] = v[i][2];
vthermal[v_dim] -= vramp;
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * vthermal[0]*vthermal[0];
t[1] += massone * vthermal[1]*vthermal[1];
t[2] += massone * vthermal[2]*vthermal[2];

12
src/compute_temp_region.cpp
View File

@ -89,18 +89,18 @@ double ComputeTempRegion::compute_scalar()
int count = 0;
double t = 0.0;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
count++;
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
count++;
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
}
}
@ -136,8 +136,8 @@ void ComputeTempRegion::compute_vector()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) {
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
t[0] += massone * v[i][0]*v[i][0];
t[1] += massone * v[i][1]*v[i][1];
t[2] += massone * v[i][2]*v[i][2];

56
src/compute_temp_sphere.cpp
View File

@ -152,7 +152,14 @@ double ComputeTempSphere::compute_scalar()
double t = 0.0;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
@ -160,13 +167,6 @@ double ComputeTempSphere::compute_scalar()
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * inertia[type[i]];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) * INERTIA*radius[i]*radius[i]*rmass[i];
}
}
if (tempbias) tbias->restore_bias_all();
@ -202,26 +202,7 @@ void ComputeTempSphere::compute_vector()
double massone,inertiaone,t[6];
for (i = 0; i < 6; i++) t[i] = 0.0;
if (mass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massone = mass[type[i]];
t[0] += massone * v[i][0]*v[i][0];
t[1] += massone * v[i][1]*v[i][1];
t[2] += massone * v[i][2]*v[i][2];
t[3] += massone * v[i][0]*v[i][1];
t[4] += massone * v[i][0]*v[i][2];
t[5] += massone * v[i][1]*v[i][2];
inertiaone = inertia[type[i]];
t[0] += inertiaone * omega[i][0]*omega[i][0];
t[1] += inertiaone * omega[i][1]*omega[i][1];
t[2] += inertiaone * omega[i][2]*omega[i][2];
t[3] += inertiaone * omega[i][0]*omega[i][1];
t[4] += inertiaone * omega[i][0]*omega[i][2];
t[5] += inertiaone * omega[i][1]*omega[i][2];
}
} else {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massone = rmass[i];
@ -240,6 +221,25 @@ void ComputeTempSphere::compute_vector()
t[4] += inertiaone * omega[i][0]*omega[i][2];
t[5] += inertiaone * omega[i][1]*omega[i][2];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massone = mass[type[i]];
t[0] += massone * v[i][0]*v[i][0];
t[1] += massone * v[i][1]*v[i][1];
t[2] += massone * v[i][2]*v[i][2];
t[3] += massone * v[i][0]*v[i][1];
t[4] += massone * v[i][0]*v[i][2];
t[5] += massone * v[i][1]*v[i][2];
inertiaone = inertia[type[i]];
t[0] += inertiaone * omega[i][0]*omega[i][0];
t[1] += inertiaone * omega[i][1]*omega[i][1];
t[2] += inertiaone * omega[i][2]*omega[i][2];
t[3] += inertiaone * omega[i][0]*omega[i][1];
t[4] += inertiaone * omega[i][0]*omega[i][2];
t[5] += inertiaone * omega[i][1]*omega[i][2];
}
}
if (tempbias) tbias->restore_bias_all();

14
src/dump_custom.cpp
View File

@ -343,15 +343,15 @@ int DumpCustom::count()
ptr = dchoose;
nstride = 1;
} else if (thresh_array[ithresh] == MASS) {
if (atom->mass) {
if (atom->rmass) {
ptr = atom->rmass;
nstride = 1;
} else {
double *mass = atom->mass;
int *type = atom->type;
for (i = 0; i < nlocal; i++) dchoose[i] = mass[type[i]];
ptr = dchoose;
nstride = 1;
} else {
ptr = atom->rmass;
nstride = 1;
}
} else if (thresh_array[ithresh] == X) {
@ -1450,16 +1450,16 @@ void DumpCustom::pack_mass(int n)
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = mass[type[i]];
buf[n] = rmass[i];
n += size_one;
}
} else {
for (int i = 0; i < nlocal; i++)
if (choose[i]) {
buf[n] = rmass[i];
buf[n] = mass[type[i]];
n += size_one;
}
}

4
src/fix_ave_spatial.cpp
View File

@ -615,8 +615,8 @@ void FixAveSpatial::end_of_step()
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (mass) values_one[layer[i]][m] += mass[type[i]];
else values_one[layer[i]][m] += rmass[i];
if (rmass) values_one[layer[i]][m] += rmass[i];
else values_one[layer[i]][m] += mass[type[i]];
// COMPUTE adds its scalar or vector component to values
// invoke compute if not previously invoked

4
src/fix_dt_reset.cpp
View File

@ -152,8 +152,8 @@ void FixDtReset::end_of_step()
if (mask[i] & groupbit) {
vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
fsq = f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
if (mass) ms = mass[type[i]];
else ms = rmass[i];
if (rmass) ms = rmass[i];
else ms = mass[type[i]];
asq = fsq/ms/ms;
bound[0] = MAX(bound[0],vsq);
bound[1] = MAX(bound[1],asq);

11
src/fix_gravity.cpp
View File

@ -152,10 +152,10 @@ void FixGravity::post_force(int vflag)
int nlocal = atom->nlocal;
double massone;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massone = mass[type[i]];
massone = rmass[i];
f[i][0] += massone*xacc;
f[i][1] += massone*yacc;
f[i][2] += massone*zacc;
@ -163,9 +163,10 @@ void FixGravity::post_force(int vflag)
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
f[i][0] += rmass[i]*xacc;
f[i][1] += rmass[i]*yacc;
f[i][2] += rmass[i]*zacc;
massone = mass[type[i]];
f[i][0] += massone*xacc;
f[i][1] += massone*yacc;
f[i][2] += massone*zacc;
}
}
}

2
src/fix_line_force.cpp
View File

@ -105,5 +105,3 @@ void FixLineForce::min_post_force(int vflag)
{
post_force(vflag);
}

12
src/fix_nve.cpp
View File

@ -74,10 +74,10 @@ void FixNVE::initial_integrate(int vflag)
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
@ -90,7 +90,7 @@ void FixNVE::initial_integrate(int vflag)
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
@ -117,10 +117,10 @@ void FixNVE::final_integrate()
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
@ -130,7 +130,7 @@ void FixNVE::final_integrate()
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];

12
src/fix_nve_limit.cpp
View File

@ -84,10 +84,10 @@ void FixNVELimit::initial_integrate(int vflag)
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
@ -110,7 +110,7 @@ void FixNVELimit::initial_integrate(int vflag)
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
@ -147,10 +147,10 @@ void FixNVELimit::final_integrate()
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
@ -169,7 +169,7 @@ void FixNVELimit::final_integrate()
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];

76
src/fix_nve_sphere.cpp
View File

@ -92,8 +92,8 @@ void FixNVESphere::init()
if (atom->mass && !atom->shape)
error->all("Fix nve/sphere requires atom attribute shape");
if (!atom->mass && (!atom->radius_flag || !atom->rmass_flag))
error->all("Fix nve/sphere requires atom attributes radius, rmass");
if (atom->rmass && !atom->radius_flag)
error->all("Fix nve/sphere requires atom attribute radius");
if (atom->mass) {
double **shape = atom->shape;
@ -137,7 +137,25 @@ void FixNVESphere::initial_integrate(int vflag)
// update v,x,omega for all particles
// d_omega/dt = torque / inertia
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
itype = type[i];
@ -155,24 +173,6 @@ void FixNVESphere::initial_integrate(int vflag)
omega[i][2] += dtirotate * torque[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
}
}
}
// update mu for dipoles
@ -228,7 +228,23 @@ void FixNVESphere::final_integrate()
dttype[i] = dtfrotate / (shape[i][0]*shape[i][0]*mass[i]);
}
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
}
}
} else {
dtfrotate = dtf / INERTIA;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
itype = type[i];
@ -243,21 +259,5 @@ void FixNVESphere::final_integrate()
omega[i][2] += dtirotate * torque[i][2];
}
}
} else {
dtfrotate = dtf / INERTIA;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
dtirotate = dtfrotate / (radius[i]*radius[i]*rmass[i]);
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
}
}
}
}

213
src/fix_shake.cpp
View File

@ -121,10 +121,10 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
type_flag[i] = 1;
} else if (mode == 'm') {
double rmass = atof(arg[next]);
if (rmass == 0.0) error->all("Invalid atom mass for fix shake");
double massone = atof(arg[next]);
if (massone == 0.0) error->all("Invalid atom mass for fix shake");
if (nmass == atom->ntypes) error->all("Too many masses for fix shake");
mass_list[nmass++] = rmass;
mass_list[nmass++] = massone;
} else error->all("Illegal fix shake command");
next++;
@ -428,6 +428,7 @@ void FixShake::pre_neighbor()
v = atom->v;
f = atom->f;
mass = atom->mass;
rmass = atom->rmass;
type = atom->type;
nlocal = atom->nlocal;
@ -562,7 +563,7 @@ void FixShake::find_clusters()
{
int i,j,m,n;
int flag,flag_all,messtag,loop,nbuf,nbufmax,size;
double imass,jmass,rmass;
double massone;
int *buf,*bufcopy;
MPI_Request request;
MPI_Status status;
@ -575,6 +576,7 @@ void FixShake::find_clusters()
int *type = atom->type;
int *mask = atom->mask;
double *mass = atom->mass;
double *rmass = atom->rmass;
int **bond_type = atom->bond_type;
int **angle_type = atom->angle_type;
int **nspecial = atom->nspecial;
@ -596,6 +598,7 @@ void FixShake::find_clusters()
// partner_tag[i][] = global IDs of each partner
// partner_mask[i][] = mask of each partner
// partner_type[i][] = type of each partner
// partner_massflag[i][] = 1 if partner meets mass criterion, 0 if not
// partner_bondtype[i][] = type of bond attached to each partner
// partner_shake[i][] = 1 if SHAKE bonded to partner, 0 if not
// partner_nshake[i][] = nshake value for each partner
@ -614,6 +617,8 @@ void FixShake::find_clusters()
memory->create_2d_int_array(nlocal,max,"shake:partner_mask");
int **partner_type =
memory->create_2d_int_array(nlocal,max,"shake:partner_type");
int **partner_massflag =
memory->create_2d_int_array(nlocal,max,"shake:partner_massflag");
int **partner_bondtype =
memory->create_2d_int_array(nlocal,max,"shake:partner_bondtype");
int **partner_shake =
@ -631,33 +636,42 @@ void FixShake::find_clusters()
}
// -----------------------------------------------------
// set partner_mask, partner_type, partner_bondtype for bonded partners
// set partner_mask, partner_type, partner_massflag, partner_bondtype
// for bonded partners
// requires communication for off-proc partners
// -----------------------------------------------------
// fill in mask, type, bondtype if own bond partner
// info to store in buf for each off-proc bond =
// 2 atoms IDs in bond, space for mask, type, bondtype
// fill in mask, type, massflag, bondtype if own bond partner
// info to store in buf for each off-proc bond = nper = 6
// 2 atoms IDs in bond, space for mask, type, massflag, bondtype
// nbufmax = largest buffer needed to hold info from any proc
int nper = 6;
nbuf = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < npartner[i]; j++) {
partner_mask[i][j] = 0;
partner_type[i][j] = 0;
partner_massflag[i][j] = 0;
partner_bondtype[i][j] = 0;
m = atom->map(partner_tag[i][j]);
if (m >= 0 && m < nlocal) {
partner_mask[i][j] = mask[m];
partner_type[i][j] = type[m];
if (nmass) {
if (rmass) massone = rmass[m];
else massone = mass[type[m]];
partner_massflag[i][j] = masscheck(massone);
}
n = bondfind(i,tag[i],partner_tag[i][j]);
if (n >= 0) partner_bondtype[i][j] = bond_type[i][n];
else {
n = bondfind(m,tag[i],partner_tag[i][j]);
if (n >= 0) partner_bondtype[i][j] = bond_type[m][n];
}
} else nbuf += 5;
} else nbuf += nper;
}
}
@ -677,18 +691,20 @@ void FixShake::find_clusters()
buf[size+1] = partner_tag[i][j];
buf[size+2] = 0;
buf[size+3] = 0;
buf[size+4] = 0;
n = bondfind(i,tag[i],partner_tag[i][j]);
if (n >= 0) buf[size+4] = bond_type[i][n];
else buf[size+4] = 0;
size += 5;
if (n >= 0) buf[size+5] = bond_type[i][n];
else buf[size+5] = 0;
size += nper;
}
}
}
// cycle buffer around ring of procs back to self
// when receive buffer, scan bond partner IDs for atoms I own
// if I own partner, fill in mask and type, search for bond with 1st atom
// and fill in bondtype
// if I own partner:
// fill in mask and type and massflag
// search for bond with 1st atom and fill in bondtype
messtag = 1;
for (loop = 0; loop < nprocs; loop++) {
@ -698,12 +714,17 @@ void FixShake::find_clusters()
if (m >= 0 && m < nlocal) {
buf[i+2] = mask[m];
buf[i+3] = type[m];
if (buf[i+4] == 0) {
if (nmass) {
if (rmass) massone = rmass[m];
else massone = mass[type[m]];
buf[i+4] = masscheck(massone);
}
if (buf[i+5] == 0) {
n = bondfind(m,buf[i],buf[i+1]);
if (n >= 0) buf[i+4] = bond_type[m][n];
if (n >= 0) buf[i+5] = bond_type[m][n];
}
}
i += 5;
i += nper;
}
if (me != next) {
MPI_Irecv(bufcopy,nbufmax,MPI_INT,prev,messtag,world,&request);
@ -723,8 +744,9 @@ void FixShake::find_clusters()
if (buf[m+1] == partner_tag[i][j]) break;
partner_mask[i][j] = buf[m+2];
partner_type[i][j] = buf[m+3];
partner_bondtype[i][j] = buf[m+4];
m += 5;
partner_massflag[i][j] = buf[m+4];
partner_bondtype[i][j] = buf[m+5];
m += nper;
}
delete [] buf;
@ -780,15 +802,17 @@ void FixShake::find_clusters()
continue;
}
if (nmass) {
imass = mass[type[i]];
jmass = mass[partner_type[i][j]];
for (m = 0; m < nmass; m++) {
rmass = mass_list[m];
if (fabs(rmass-imass) <= MASSDELTA ||
fabs(rmass-jmass) <= MASSDELTA) {
if (partner_massflag[i][j]) {
partner_shake[i][j] = 1;
nshake[i]++;
continue;
} else {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
if (masscheck(massone)) {
partner_shake[i][j] = 1;
nshake[i]++;
break;
continue;
}
}
}
@ -1049,6 +1073,7 @@ void FixShake::find_clusters()
memory->destroy_2d_int_array(partner_tag);
memory->destroy_2d_int_array(partner_mask);
memory->destroy_2d_int_array(partner_type);
memory->destroy_2d_int_array(partner_massflag);
memory->destroy_2d_int_array(partner_bondtype);
memory->destroy_2d_int_array(partner_shake);
memory->destroy_2d_int_array(partner_nshake);
@ -1124,6 +1149,18 @@ void FixShake::find_clusters()
}
}
/* ----------------------------------------------------------------------
check if massone is within MASSDELTA of any mass in mass_list
return 1 if yes, 0 if not
------------------------------------------------------------------------- */
int FixShake::masscheck(double massone)
{
for (int i = 0; i < nmass; i++)
if (fabs(mass_list[i]-massone) <= MASSDELTA) return 1;
return 0;
}
/* ----------------------------------------------------------------------
update the unconstrained position of each atom
only for SHAKE clusters, else set to 0.0
@ -1134,13 +1171,24 @@ void FixShake::unconstrained_update()
{
double dtfmsq;
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
dtfmsq = dtfsq / mass[type[i]];
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
dtfmsq = dtfsq / rmass[i];
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
}
} else {
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
dtfmsq = dtfsq / mass[type[i]];
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
}
}
}
@ -1150,6 +1198,7 @@ void FixShake::shake2(int m)
{
int nlist,list[2];
double v[6];
double invmass0,invmass1;
// local atom IDs and constraint distances
@ -1180,8 +1229,13 @@ void FixShake::shake2(int m)
// a,b,c = coeffs in quadratic equation for lamda
double invmass0 = 1.0/mass[type[i0]];
double invmass1 = 1.0/mass[type[i1]];
if (rmass) {
invmass0 = 1.0/rmass[i0];
invmass1 = 1.0/rmass[i1];
} else {
invmass0 = 1.0/mass[type[i0]];
invmass1 = 1.0/mass[type[i1]];
}
double a = (invmass0+invmass1)*(invmass0+invmass1) * r01sq;
double b = 2.0 * (invmass0+invmass1) *
@ -1243,6 +1297,7 @@ void FixShake::shake3(int m)
{
int nlist,list[3];
double v[6];
double invmass0,invmass1,invmass2;
// local atom IDs and constraint distances
@ -1289,9 +1344,15 @@ void FixShake::shake3(int m)
// matrix coeffs and rhs for lamda equations
double invmass0 = 1.0/mass[type[i0]];
double invmass1 = 1.0/mass[type[i1]];
double invmass2 = 1.0/mass[type[i2]];
if (rmass) {
invmass0 = 1.0/rmass[i0];
invmass1 = 1.0/rmass[i1];
invmass2 = 1.0/rmass[i2];
} else {
invmass0 = 1.0/mass[type[i0]];
invmass1 = 1.0/mass[type[i1]];
invmass2 = 1.0/mass[type[i2]];
}
double a11 = 2.0 * (invmass0+invmass1) *
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
@ -1401,6 +1462,7 @@ void FixShake::shake4(int m)
{
int nlist,list[4];
double v[6];
double invmass0,invmass1,invmass2,invmass3;
// local atom IDs and constraint distances
@ -1463,10 +1525,17 @@ void FixShake::shake4(int m)
// matrix coeffs and rhs for lamda equations
double invmass0 = 1.0/mass[type[i0]];
double invmass1 = 1.0/mass[type[i1]];
double invmass2 = 1.0/mass[type[i2]];
double invmass3 = 1.0/mass[type[i3]];
if (rmass) {
invmass0 = 1.0/rmass[i0];
invmass1 = 1.0/rmass[i1];
invmass2 = 1.0/rmass[i2];
invmass3 = 1.0/rmass[i3];
} else {
invmass0 = 1.0/mass[type[i0]];
invmass1 = 1.0/mass[type[i1]];
invmass2 = 1.0/mass[type[i2]];
invmass3 = 1.0/mass[type[i3]];
}
double a11 = 2.0 * (invmass0+invmass1) *
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
@ -1636,6 +1705,7 @@ void FixShake::shake3angle(int m)
{
int nlist,list[3];
double v[6];
double invmass0,invmass1,invmass2;
// local atom IDs and constraint distances
@ -1697,9 +1767,15 @@ void FixShake::shake3angle(int m)
// matrix coeffs and rhs for lamda equations
double invmass0 = 1.0/mass[type[i0]];
double invmass1 = 1.0/mass[type[i1]];
double invmass2 = 1.0/mass[type[i2]];
if (rmass) {
invmass0 = 1.0/rmass[i0];
invmass1 = 1.0/rmass[i1];
invmass2 = 1.0/rmass[i2];
} else {
invmass0 = 1.0/mass[type[i0]];
invmass1 = 1.0/mass[type[i1]];
invmass2 = 1.0/mass[type[i2]];
}
double a11 = 2.0 * (invmass0+invmass1) *
(s01[0]*r01[0] + s01[1]*r01[1] + s01[2]*r01[2]);
@ -2204,20 +2280,39 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
double invmass,dtfmsq;
int jlevel;
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
invmass = 1.0 / mass[type[i]];
dtfmsq = dtfsq * invmass;
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
for (jlevel = 0; jlevel < ilevel; jlevel++) {
dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass;
xshake[i][0] += dtfmsq*f_level[i][jlevel][0];
xshake[i][1] += dtfmsq*f_level[i][jlevel][1];
xshake[i][2] += dtfmsq*f_level[i][jlevel][2];
}
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
invmass = 1.0 / rmass[i];
dtfmsq = dtfsq * invmass;
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
for (jlevel = 0; jlevel < ilevel; jlevel++) {
dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass;
xshake[i][0] += dtfmsq*f_level[i][jlevel][0];
xshake[i][1] += dtfmsq*f_level[i][jlevel][1];
xshake[i][2] += dtfmsq*f_level[i][jlevel][2];
}
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
}
} else {
for (int i = 0; i < nlocal; i++) {
if (shake_flag[i]) {
invmass = 1.0 / mass[type[i]];
dtfmsq = dtfsq * invmass;
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0];
xshake[i][1] = x[i][1] + dtv*v[i][1] + dtfmsq*f[i][1];
xshake[i][2] = x[i][2] + dtv*v[i][2] + dtfmsq*f[i][2];
for (jlevel = 0; jlevel < ilevel; jlevel++) {
dtfmsq = dtf_innerhalf * step_respa[jlevel] * invmass;
xshake[i][0] += dtfmsq*f_level[i][jlevel][0];
xshake[i][1] += dtfmsq*f_level[i][jlevel][1];
xshake[i][2] += dtfmsq*f_level[i][jlevel][2];
}
} else xshake[i][2] = xshake[i][1] = xshake[i][0] = 0.0;
}
}
// communicate results if necessary

5
src/fix_shake.h
View File

@ -52,7 +52,7 @@ class FixShake : public Fix {
int *bond_flag,*angle_flag; // bond/angle types to constrain
int *type_flag; // constrain bonds to these types
double *mass_list; // constrain bonds to these masses
int nmass;
int nmass; // # of masses in mass_list
double *bond_distance,*angle_distance; // constraint distances
@ -62,7 +62,7 @@ class FixShake : public Fix {
double *step_respa;
double **x,**v,**f; // local ptrs to atom class quantities
double *mass;
double *mass,*rmass;
int *type;
int nlocal;
// atom-based arrays
@ -93,6 +93,7 @@ class FixShake : public Fix {
double *a_ave_all,*a_max_all,*a_min_all;
void find_clusters();
int masscheck(double);
void unconstrained_update();
void shake2(int);
void shake3(int);

8
src/fix_thermal_conductivity.cpp
View File

@ -182,8 +182,8 @@ void FixThermalConductivity::end_of_step()
if (coord >= slablo_lo && coord < slablo_hi) {
ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (mass) ke *= 0.5*mass[type[i]];
else ke *= 0.5*rmass[i];
if (rmass) ke *= 0.5*rmass[i];
else ke *= 0.5*mass[type[i]];
if (nlo < nswap || ke > ke_lo[nswap-1]) {
for (insert = nlo-1; insert >= 0; insert--)
if (ke < ke_lo[insert]) break;
@ -200,8 +200,8 @@ void FixThermalConductivity::end_of_step()
if (coord >= slabhi_lo && coord < slabhi_hi) {
ke = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (mass) ke *= 0.5*mass[type[i]];
else ke *= 0.5*rmass[i];
if (rmass) ke *= 0.5*rmass[i];
else ke *= 0.5*mass[type[i]];
if (nhi < nswap || ke < ke_hi[nswap-1]) {
for (insert = nhi-1; insert >= 0; insert--)
if (ke > ke_hi[insert]) break;

19
src/fix_viscosity.cpp
View File

@ -176,8 +176,6 @@ void FixViscosity::end_of_step()
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -232,6 +230,9 @@ void FixViscosity::end_of_step()
// exchange momenta between the 2 particles
// if I own both particles just swap, else point2point comm of vel,mass
double *mass = atom->mass;
double *rmass = atom->rmass;
int ipos,ineg;
double sbuf[2],rbuf[2];
@ -254,10 +255,10 @@ void FixViscosity::end_of_step()
ipos = pos_index[ipositive++];
ineg = neg_index[inegative++];
rbuf[0] = v[ipos][vdim];
if (mass) rbuf[1] = mass[type[ipos]];
else rbuf[1] = rmass[ipos];
if (rmass) rbuf[1] = rmass[ipos];
else rbuf[1] = mass[type[ipos]];
sbuf[0] = v[ineg][vdim];
if (mass) sbuf[1] = mass[type[ineg]];
if (rmass) sbuf[1] = mass[type[ineg]];
else sbuf[1] = rmass[ineg];
v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1];
v[ipos][vdim] = sbuf[0] * sbuf[1]/rbuf[1];
@ -266,8 +267,8 @@ void FixViscosity::end_of_step()
} else if (me == all[0].proc) {
ipos = pos_index[ipositive++];
sbuf[0] = v[ipos][vdim];
if (mass) sbuf[1] = mass[type[ipos]];
else sbuf[1] = rmass[ipos];
if (rmass) sbuf[1] = rmass[ipos];
else sbuf[1] = mass[type[ipos]];
MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0,
rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status);
v[ipos][vdim] = rbuf[0] * rbuf[1]/sbuf[1];
@ -276,8 +277,8 @@ void FixViscosity::end_of_step()
} else if (me == all[1].proc) {
ineg = neg_index[inegative++];
sbuf[0] = v[ineg][vdim];
if (mass) sbuf[1] = mass[type[ineg]];
else sbuf[1] = rmass[ineg];
if (rmass) sbuf[1] = rmass[ineg];
else sbuf[1] = mass[type[ineg]];
MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0,
rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status);
v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1];

44
src/group.cpp
View File

@ -446,12 +446,12 @@ double Group::mass(int igroup)
double one = 0.0;
if (mass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) one += mass[type[i]];
} else {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) one += rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) one += mass[type[i]];
}
double all;
@ -550,13 +550,13 @@ void Group::xcm(int igroup, double masstotal, double *cm)
int xbox,ybox,zbox;
double massone;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
massone = mass[type[i]];
massone = rmass[i];
cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone;
@ -567,7 +567,7 @@ void Group::xcm(int igroup, double masstotal, double *cm)
xbox = (image[i] & 1023) - 512;
ybox = (image[i] >> 10 & 1023) - 512;
zbox = (image[i] >> 20) - 512;
massone = rmass[i];
massone = mass[type[i]];
cmone[0] += (x[i][0] + xbox*xprd) * massone;
cmone[1] += (x[i][1] + ybox*yprd) * massone;
cmone[2] += (x[i][2] + zbox*zprd) * massone;
@ -599,10 +599,10 @@ void Group::vcm(int igroup, double masstotal, double *cm)
double p[3],massone;
p[0] = p[1] = p[2] = 0.0;
if (mass) {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massone = mass[type[i]];
massone = rmass[i];
p[0] += v[i][0]*massone;
p[1] += v[i][1]*massone;
p[2] += v[i][2]*massone;
@ -610,7 +610,7 @@ void Group::vcm(int igroup, double masstotal, double *cm)
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massone = rmass[i];
massone = mass[type[i]];
p[0] += v[i][0]*massone;
p[1] += v[i][1]*massone;
p[2] += v[i][2]*massone;
@ -665,16 +665,16 @@ double Group::ke(int igroup)
double one = 0.0;
if (mass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
} else {
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
rmass[i];
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
one += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
}
double all;
@ -715,8 +715,8 @@ double Group::gyration(int igroup, double masstotal, double *cm)
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
rg += (dx*dx + dy*dy + dz*dz) * massone;
}
double rg_all;
@ -759,8 +759,8 @@ void Group::angmom(int igroup, double *cm, double *lmom)
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
p[0] += massone * (dy*v[i][2] - dz*v[i][1]);
p[1] += massone * (dz*v[i][0] - dx*v[i][2]);
p[2] += massone * (dx*v[i][1] - dy*v[i][0]);
@ -806,8 +806,8 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3])
dx = (x[i][0] + xbox*xprd) - cm[0];
dy = (x[i][1] + ybox*yprd) - cm[1];
dz = (x[i][2] + zbox*zprd) - cm[2];
if (mass) massone = mass[type[i]];
else massone = rmass[i];
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
ione[0][0] += massone * (dy*dy + dz*dz);
ione[1][1] += massone * (dx*dx + dz*dz);
ione[2][2] += massone * (dx*dx + dy*dy);

8
src/pair_lj_cut.cpp
View File

@ -714,11 +714,3 @@ double PairLJCut::single(int i, int j, int itype, int jtype, double rsq,
offset[itype][jtype];
return factor_lj*philj;
}
/* ---------------------------------------------------------------------- */
void *PairLJCut::extract(char *str)
{
if (strcmp(str,"diameter") == 0) return (void *) sigma;
return NULL;
}

1
src/pair_lj_cut.h
View File

@ -37,7 +37,6 @@ class PairLJCut : public Pair {
void compute_inner();
void compute_middle();
void compute_outer(int, int);
void *extract(char *);
protected:
double cut_global;

12
src/velocity.cpp
View File

@ -225,8 +225,8 @@ void Velocity::create(int narg, char **arg)
m = atom->map(i);
if (m >= 0 && m < nlocal) {
if (mask[m] & groupbit) {
if (mass) factor = 1.0/sqrt(mass[type[m]]);
else factor = 1.0/sqrt(rmass[m]);
if (rmass) factor = 1.0/sqrt(rmass[m]);
factor = 1.0/sqrt(mass[type[m]]);
v[m][0] = vx * factor;
v[m][1] = vy * factor;
if (dimension == 3) v[m][2] = vz * factor;
@ -257,8 +257,8 @@ void Velocity::create(int narg, char **arg)
vy = random->gaussian();
vz = random->gaussian();
}
if (mass) factor = 1.0/sqrt(mass[type[i]]);
else factor = 1.0/sqrt(rmass[i]);
if (rmass) factor = 1.0/sqrt(rmass[i]);
factor = 1.0/sqrt(mass[type[i]]);
v[i][0] = vx * factor;
v[i][1] = vy * factor;
if (dimension == 3) v[i][2] = vz * factor;
@ -283,8 +283,8 @@ void Velocity::create(int narg, char **arg)
vz = random->gaussian();
}
if (mass) factor = 1.0/sqrt(mass[type[i]]);
else factor = 1.0/sqrt(rmass[i]);
if (rmass) factor = 1.0/sqrt(rmass[i]);
factor = 1.0/sqrt(mass[type[i]]);
v[i][0] = vx * factor;
v[i][1] = vy * factor;
if (dimension == 3) v[i][2] = vz * factor;