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

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

Browse Source
This commit is contained in:
sjplimp
2009-01-05 22:26:08 +00:00
parent e1f72f6d3a
commit 0cc1d090f4
90 changed files with 1129 additions and 779 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/ASPHERE/compute_erotate_asphere.cpp
View File

@ -15,14 +15,13 @@
#include "compute_erotate_asphere.h"
#include "math_extra.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
/* ---------------------------------------------------------------------- */
ComputeERotateAsphere::ComputeERotateAsphere(LAMMPS *lmp, int narg, char **arg) :
@ -63,7 +62,7 @@ void ComputeERotateAsphere::init()
double ComputeERotateAsphere::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **quat = atom->quat;
double **angmom = atom->angmom;

13
src/ASPHERE/compute_temp_asphere.cpp
View File

@ -19,6 +19,7 @@
#include "compute_temp_asphere.h"
#include "math_extra.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
@ -29,9 +30,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
@ -150,11 +148,10 @@ void ComputeTempAsphere::dof_compute()
double ComputeTempAsphere::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
if (tempbias) {
if (!(tbias->invoked & INVOKED_SCALAR))
double tmp = tbias->compute_scalar();
if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar();
tbias->remove_bias_all();
}
@ -213,10 +210,10 @@ void ComputeTempAsphere::compute_vector()
{
int i;
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
if (tempbias) {
if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector();
if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector();
tbias->remove_bias_all();
}

2
src/GRANULAR/fix_pour.cpp
View File

@ -42,6 +42,8 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 6) error->all("Illegal fix pour command");
time_depend = 1;
if (!atom->radius_flag || !atom->rmass_flag)
error->all("Fix pour requires atom attributes radius, rmass");

227
src/GRANULAR/fix_wall_gran.cpp
View File

@ -45,6 +45,8 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 4) error->all("Illegal fix wall/gran command");
time_depend = 1;
if (!atom->radius_flag || !atom->omega_flag || !atom->torque_flag)
error->all("Fix wall/gran requires atom attributes radius, omega, torque");
@ -199,10 +201,8 @@ void FixWallGran::init()
// same initialization as in pair_gran_history::init_style()
xkkt = xkk * 2.0/7.0;
dt = update->dt;
double gammas = 0.5*gamman;
gamman_dl = gamman/dt;
gammas_dl = gammas/dt;
dt = update->dt;
// set pairstyle from granular pair style
@ -324,9 +324,11 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
{
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,ft,fs1,fs2,fs3,ccelx,ccely,ccelz,tor1,tor2,tor3,rinv;
double fn,fs,ft,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3,rinv,rsqinv;
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
@ -334,16 +336,12 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
vr1 *= dt;
vr2 *= dt;
vr3 *= dt;
// normal component
vnnr = vr1*dx + vr2*dy + vr3*dz;
vn1 = dx*vnnr / rsq;
vn2 = dy*vnnr / rsq;
vn3 = dz*vnnr / rsq;
vn1 = dx*vnnr * rsqinv;
vn2 = dy*vnnr * rsqinv;
vn3 = dz*vnnr * rsqinv;
// tangential component
@ -353,20 +351,15 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
// relative rotational velocity
wr1 = radius*omega[0];
wr2 = radius*omega[1];
wr3 = radius*omega[2];
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
wr1 *= dt/r;
wr2 *= dt/r;
wr3 *= dt/r;
// normal damping term
// this definition of DAMP includes the extra 1/r term
// normal forces = Hookian contact + normal velocity damping
xmeff = mass;
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radius-r)/r - damp;
damp = xmeff*gamman*vnnr*rsqinv;
ccel = xkk*(radius-r)*rinv - damp;
// relative velocities
@ -379,31 +372,26 @@ void FixWallGran::no_history(double rsq, double dx, double dy, double dz,
// force normalization
fn = xmu * fabs(ccel*r);
fs = xmeff*gammas_dl*vrel;
fs = xmeff*gammas*vrel;
if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
else ft = 0.0;
// shear friction forces
// tangential force due to tangential velocity damping
fs1 = -ft*vtr1;
fs2 = -ft*vtr2;
fs3 = -ft*vtr3;
// force components
// forces & torques
ccelx = dx*ccel + fs1;
ccely = dy*ccel + fs2;
ccelz = dz*ccel + fs3;
fx = dx*ccel + fs1;
fy = dy*ccel + fs2;
fz = dz*ccel + fs3;
// forces
f[0] += fx;
f[1] += fy;
f[2] += fz;
f[0] += ccelx;
f[1] += ccely;
f[2] += ccelz;
// torques
rinv = 1/r;
tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1);
@ -421,10 +409,12 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,
{
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,ccelx,ccely,ccelz,tor1,tor2,tor3;
double shrmag,rsht,rinv;
double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
double shrmag,rsht,rinv,rsqinv;
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
@ -432,16 +422,12 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
vr1 *= dt;
vr2 *= dt;
vr3 *= dt;
// normal component
vnnr = vr1*dx + vr2*dy + vr3*dz;
vn1 = dx*vnnr / rsq;
vn2 = dy*vnnr / rsq;
vn3 = dz*vnnr / rsq;
vn1 = dx*vnnr * rsqinv;
vn2 = dy*vnnr * rsqinv;
vn3 = dz*vnnr * rsqinv;
// tangential component
@ -451,20 +437,15 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,
// relative rotational velocity
wr1 = radius*omega[0];
wr2 = radius*omega[1];
wr3 = radius*omega[2];
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
wr1 *= dt/r;
wr2 *= dt/r;
wr3 *= dt/r;
// normal damping term
// this definition of DAMP includes the extra 1/r term
// normal forces = Hookian contact + normal velocity damping
xmeff = mass;
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radius-r)/r - damp;
damp = xmeff*gamman*vnnr*rsqinv;
ccel = xkk*(radius-r)*rinv - damp;
// relative velocities
@ -476,60 +457,54 @@ void FixWallGran::history(double rsq, double dx, double dy, double dz,
// shear history effects
shear[0] += vtr1;
shear[1] += vtr2;
shear[2] += vtr3;
shear[0] += vtr1*dt;
shear[1] += vtr2*dt;
shear[2] += vtr3*dt;
shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]);
// rotate shear displacements correctly
// rotate shear displacements
rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz;
rsht = rsht/rsq;
rsht = rsht*rsqinv;
shear[0] -= rsht*dx;
shear[1] -= rsht*dy;
shear[2] -= rsht*dz;
// tangential forces
// tangential forces = shear + tangential velocity damping
fs1 = - (xkkt*shear[0] + xmeff*gammas_dl*vtr1);
fs2 = - (xkkt*shear[1] + xmeff*gammas_dl*vtr2);
fs3 = - (xkkt*shear[2] + xmeff*gammas_dl*vtr3);
fs1 = - (xkkt*shear[0] + xmeff*gammas*vtr1);
fs2 = - (xkkt*shear[1] + xmeff*gammas*vtr2);
fs3 = - (xkkt*shear[2] + xmeff*gammas*vtr3);
// force normalization
// rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
fn = xmu * fabs(ccel*r);
// shrmag is magnitude of shearwall
// rescale frictional displacements and forces if needed
if (fs > fn) {
if (shrmag != 0.0) {
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) -
xmeff*gammas_dl*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) -
xmeff*gammas_dl*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) -
xmeff*gammas_dl*vtr3/xkkt;
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) -
xmeff*gammas*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) -
xmeff*gammas*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) -
xmeff*gammas*vtr3/xkkt;
fs1 = fs1 * fn / fs ;
fs2 = fs2 * fn / fs;
fs3 = fs3 * fn / fs;
} else fs1 = fs2 = fs3 = 0.0;
}
ccelx = dx*ccel + fs1;
ccely = dy*ccel + fs2;
ccelz = dz*ccel + fs3;
// forces & torques
// forces
fx = dx*ccel + fs1;
fy = dy*ccel + fs2;
fz = dz*ccel + fs3;
f[0] += ccelx;
f[1] += ccely;
f[2] += ccelz;
f[0] += fx;
f[1] += fy;
f[2] += fz;
// torques
rinv = 1/r;
tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1);
@ -547,10 +522,12 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
{
double r,vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3,xmeff,damp,ccel,vtr1,vtr2,vtr3,vrel;
double fn,fs,fs1,fs2,fs3,ccelx,ccely,ccelz,tor1,tor2,tor3;
double shrmag,rsht,rhertz,rinv;
double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3;
double shrmag,rsht,rhertz,rinv,rsqinv;
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
@ -558,10 +535,6 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
vr2 = v[1] - vwall[1];
vr3 = v[2] - vwall[2];
vr1 *= dt;
vr2 *= dt;
vr3 *= dt;
// normal component
vnnr = vr1*dx + vr2*dy + vr3*dz;
@ -577,20 +550,15 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
// relative rotational velocity
wr1 = radius*omega[0];
wr2 = radius*omega[1];
wr3 = radius*omega[2];
wr1 = radius*omega[0] * rinv;
wr2 = radius*omega[1] * rinv;
wr3 = radius*omega[2] * rinv;
wr1 *= dt/r;
wr2 *= dt/r;
wr3 *= dt/r;
// normal damping term
// this definition of DAMP includes the extra 1/r term
// normal forces = Hertzian contact + normal velocity damping
xmeff = mass;
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radius-r)/r - damp;
damp = xmeff*gamman*vnnr*rsqinv;
ccel = xkk*(radius-r)*rinv - damp;
rhertz = sqrt(radius - r);
ccel = rhertz * ccel;
@ -604,60 +572,54 @@ void FixWallGran::hertzian(double rsq, double dx, double dy, double dz,
// shear history effects
shear[0] += vtr1;
shear[1] += vtr2;
shear[2] += vtr3;
shear[0] += vtr1*dt;
shear[1] += vtr2*dt;
shear[2] += vtr3*dt;
shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]);
// rotate shear displacements correctly
// rotate shear displacements
rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz;
rsht = rsht/rsq;
rsht = rsht*rsqinv;
shear[0] -= rsht*dx;
shear[1] -= rsht*dy;
shear[2] -= rsht*dz;
// tangential forces
// tangential forces = shear + tangential velocity damping
fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas_dl*vtr1);
fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas_dl*vtr2);
fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas_dl*vtr3);
fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas*vtr1);
fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas*vtr2);
fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas*vtr3);
// force normalization
// rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
fn = xmu * fabs(ccel*r);
// shrmag is magnitude of shearwall
// rescale frictional displacements and forces if needed
if (fs > fn) {
if (shrmag != 0.0) {
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) -
xmeff*gammas_dl*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) -
xmeff*gammas_dl*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) -
xmeff*gammas_dl*vtr3/xkkt;
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) -
xmeff*gammas*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) -
xmeff*gammas*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) -
xmeff*gammas*vtr3/xkkt;
fs1 = fs1 * fn / fs ;
fs2 = fs2 * fn / fs;
fs3 = fs3 * fn / fs;
} else fs1 = fs2 = fs3 = 0.0;
}
ccelx = dx*ccel + fs1;
ccely = dy*ccel + fs2;
ccelz = dz*ccel + fs3;
// forces & torques
// forces
fx = dx*ccel + fs1;
fy = dy*ccel + fs2;
fz = dz*ccel + fs3;
f[0] += ccelx;
f[1] += ccely;
f[2] += ccelz;
f[0] += fx;
f[1] += fy;
f[2] += fz;
// torques
rinv = 1/r;
tor1 = rinv * (dy*fs3 - dz*fs2);
tor2 = rinv * (dz*fs1 - dx*fs3);
tor3 = rinv * (dx*fs2 - dy*fs1);
@ -778,7 +740,4 @@ int FixWallGran::size_restart(int nlocal)
void FixWallGran::reset_dt()
{
dt = update->dt;
double gammas = 0.5*gamman;
gamman_dl = gamman/dt;
gammas_dl = gammas/dt;
}

4
src/GRANULAR/fix_wall_gran.h
View File

@ -40,9 +40,9 @@ class FixWallGran : public Fix {
private:
int wallstyle,pairstyle,wiggle,wshear,axis;
double xkk,xkkt,gamman,xmu;
double xkk,xkkt,gamman,gammas,xmu;
double lo,hi,cylradius;
double dt,gamman_dl,gammas_dl;
double dt;
double amplitude,period,omega,time_origin,vshear;
int *touch;

66
src/GRANULAR/pair_gran_hertzian.cpp
View File

@ -42,7 +42,7 @@ void PairGranHertzian::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
double radi,radj,radsum,rsq,r,rinv;
double radi,radj,radsum,rsq,r,rinv,rsqinv;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
@ -108,6 +108,8 @@ void PairGranHertzian::compute(int eflag, int vflag)
} else {
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
@ -115,16 +117,12 @@ void PairGranHertzian::compute(int eflag, int vflag)
vr2 = v[i][1] - v[j][1];
vr3 = v[i][2] - v[j][2];
vr1 *= dt;
vr2 *= dt;
vr3 *= dt;
// normal component
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr / rsq;
vn2 = dely*vnnr / rsq;
vn3 = delz*vnnr / rsq;
vn1 = delx*vnnr * rsqinv;
vn2 = dely*vnnr * rsqinv;
vn3 = delz*vnnr * rsqinv;
// tangential component
@ -134,24 +132,19 @@ void PairGranHertzian::compute(int eflag, int vflag)
// relative rotational velocity
wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
wr1 *= dt/r;
wr2 *= dt/r;
wr3 *= dt/r;
// normal damping term
// this definition of DAMP includes the extra 1/r term
// normal force = Hertzian contact + normal velocity damping
xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) xmeff = rmass[j];
if (mask[j] & freeze_group_bit) xmeff = rmass[i];
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radsum-r)/r - damp;
damp = xmeff*gamman*vnnr*rsqinv;
ccel = xkk*(radsum-r)*rinv - damp;
rhertz = sqrt(radsum - r);
ccel = rhertz * ccel;
ccel *= rhertz;
// relative velocities
@ -162,31 +155,29 @@ void PairGranHertzian::compute(int eflag, int vflag)
vrel = sqrt(vrel);
// shear history effects
// shrmag = magnitude of shear
touch[jj] = 1;
shear = &allshear[3*jj];
shear[0] += vtr1;
shear[1] += vtr2;
shear[2] += vtr3;
shear[0] += vtr1*dt;
shear[1] += vtr2*dt;
shear[2] += vtr3*dt;
shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] +
shear[2]*shear[2]);
// rotate shear displacements correctly
// rotate shear displacements
rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz;
rsht /= rsq;
rsht *= rsqinv;
shear[0] -= rsht*delx;
shear[1] -= rsht*dely;
shear[2] -= rsht*delz;
// tangential forces
// tangential forces = shear + tangential velocity damping
fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas_dl*vtr1);
fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas_dl*vtr2);
fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas_dl*vtr3);
fs1 = -rhertz * (xkkt*shear[0] + xmeff*gammas*vtr1);
fs2 = -rhertz * (xkkt*shear[1] + xmeff*gammas*vtr2);
fs3 = -rhertz * (xkkt*shear[2] + xmeff*gammas*vtr3);
// force normalization
// rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
@ -194,12 +185,12 @@ void PairGranHertzian::compute(int eflag, int vflag)
if (fs > fn) {
if (shrmag != 0.0) {
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) -
xmeff*gammas_dl*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) -
xmeff*gammas_dl*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) -
xmeff*gammas_dl*vtr3/xkkt;
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) -
xmeff*gammas*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) -
xmeff*gammas*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) -
xmeff*gammas*vtr3/xkkt;
fs1 *= fn/fs;
fs2 *= fn/fs;
fs3 *= fn/fs;
@ -219,7 +210,6 @@ void PairGranHertzian::compute(int eflag, int vflag)
f[i][1] += fy;
f[i][2] += fz;
rinv = 1/r;
tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1);

98
src/GRANULAR/pair_gran_history.cpp
View File

@ -69,7 +69,7 @@ void PairGranHistory::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
double radi,radj,radsum,rsq,r,rinv;
double radi,radj,radsum,rsq,r,rinv,rsqinv;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
@ -135,6 +135,8 @@ void PairGranHistory::compute(int eflag, int vflag)
} else {
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
@ -142,16 +144,12 @@ void PairGranHistory::compute(int eflag, int vflag)
vr2 = v[i][1] - v[j][1];
vr3 = v[i][2] - v[j][2];
vr1 *= dt;
vr2 *= dt;
vr3 *= dt;
// normal component
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr / rsq;
vn2 = dely*vnnr / rsq;
vn3 = delz*vnnr / rsq;
vn1 = delx*vnnr * rsqinv;
vn2 = dely*vnnr * rsqinv;
vn3 = delz*vnnr * rsqinv;
// tangential component
@ -161,22 +159,17 @@ void PairGranHistory::compute(int eflag, int vflag)
// relative rotational velocity
wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
wr1 *= dt/r;
wr2 *= dt/r;
wr3 *= dt/r;
// normal damping term
// this definition of DAMP includes the extra 1/r term
// normal forces = Hookian contact + normal velocity damping
xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) xmeff = rmass[j];
if (mask[j] & freeze_group_bit) xmeff = rmass[i];
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radsum-r)/r - damp;
damp = xmeff*gamman*vnnr*rsqinv;
ccel = xkk*(radsum-r)*rinv - damp;
// relative velocities
@ -187,31 +180,29 @@ void PairGranHistory::compute(int eflag, int vflag)
vrel = sqrt(vrel);
// shear history effects
// shrmag = magnitude of shear
touch[jj] = 1;
shear = &allshear[3*jj];
shear[0] += vtr1;
shear[1] += vtr2;
shear[2] += vtr3;
shear[0] += vtr1*dt;
shear[1] += vtr2*dt;
shear[2] += vtr3*dt;
shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] +
shear[2]*shear[2]);
// rotate shear displacements correctly
// rotate shear displacements
rsht = shear[0]*delx + shear[1]*dely + shear[2]*delz;
rsht /= rsq;
rsht *= rsqinv;
shear[0] -= rsht*delx;
shear[1] -= rsht*dely;
shear[2] -= rsht*delz;
// tangential forces
// tangential forces = shear + tangential velocity damping
fs1 = - (xkkt*shear[0] + xmeff*gammas_dl*vtr1);
fs2 = - (xkkt*shear[1] + xmeff*gammas_dl*vtr2);
fs3 = - (xkkt*shear[2] + xmeff*gammas_dl*vtr3);
fs1 = - (xkkt*shear[0] + xmeff*gammas*vtr1);
fs2 = - (xkkt*shear[1] + xmeff*gammas*vtr2);
fs3 = - (xkkt*shear[2] + xmeff*gammas*vtr3);
// force normalization
// rescale frictional displacements and forces if needed
fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
@ -219,12 +210,12 @@ void PairGranHistory::compute(int eflag, int vflag)
if (fs > fn) {
if (shrmag != 0.0) {
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas_dl*vtr1/xkkt) -
xmeff*gammas_dl*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas_dl*vtr2/xkkt) -
xmeff*gammas_dl*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas_dl*vtr3/xkkt) -
xmeff*gammas_dl*vtr3/xkkt;
shear[0] = (fn/fs) * (shear[0] + xmeff*gammas*vtr1/xkkt) -
xmeff*gammas*vtr1/xkkt;
shear[1] = (fn/fs) * (shear[1] + xmeff*gammas*vtr2/xkkt) -
xmeff*gammas*vtr2/xkkt;
shear[2] = (fn/fs) * (shear[2] + xmeff*gammas*vtr3/xkkt) -
xmeff*gammas*vtr3/xkkt;
fs1 *= fn/fs;
fs2 *= fn/fs;
fs3 *= fn/fs;
@ -244,7 +235,6 @@ void PairGranHistory::compute(int eflag, int vflag)
f[i][1] += fy;
f[i][2] += fz;
rinv = 1/r;
tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1);
@ -297,15 +287,6 @@ void PairGranHistory::settings(int narg, char **arg)
gamman = atof(arg[1]);
xmu = atof(arg[2]);
dampflag = atoi(arg[3]);
// granular styles do not use pair_coeff, so set setflag for everything now
if (!allocated) allocate();
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++)
setflag[i][j] = 1;
}
/* ----------------------------------------------------------------------
@ -314,7 +295,22 @@ void PairGranHistory::settings(int narg, char **arg)
void PairGranHistory::coeff(int narg, char **arg)
{
error->all("Granular pair styles do not use pair_coeff settings");
if (narg > 2) error->all("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);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all("Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
@ -344,11 +340,9 @@ void PairGranHistory::init_style()
}
xkkt = xkk * 2.0/7.0;
dt = update->dt;
double gammas = 0.5*gamman;
if (dampflag == 0) gammas = 0.0;
gamman_dl = gamman/dt;
gammas_dl = gammas/dt;
dt = update->dt;
// if shear history is stored:
// check if newton flag is valid
@ -512,8 +506,4 @@ void *PairGranHistory::extract(char *str)
void PairGranHistory::reset_dt()
{
dt = update->dt;
double gammas = 0.5*gamman;
if (dampflag == 0) gammas = 0.0;
gamman_dl = gamman/dt;
gammas_dl = gammas/dt;
}

4
src/GRANULAR/pair_gran_history.h
View File

@ -38,8 +38,8 @@ class PairGranHistory : public Pair {
protected:
double xkk,xkkt,xmu;
int dampflag;
double gamman;
double dt,gamman_dl,gammas_dl;
double gamman,gammas;
double dt;
int freeze_group_bit;
int history;
class FixShearHistory *fix_history;

36
src/GRANULAR/pair_gran_no_history.cpp
View File

@ -41,7 +41,7 @@ void PairGranNoHistory::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
double radi,radj,radsum,rsq,r,rinv;
double radi,radj,radsum,rsq,r,rinv,rsqinv;
double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
double wr1,wr2,wr3;
double vtr1,vtr2,vtr3,vrel;
@ -91,6 +91,8 @@ void PairGranNoHistory::compute(int eflag, int vflag)
if (rsq < radsum*radsum) {
r = sqrt(rsq);
rinv = 1.0/r;
rsqinv = 1.0/rsq;
// relative translational velocity
@ -98,16 +100,12 @@ void PairGranNoHistory::compute(int eflag, int vflag)
vr2 = v[i][1] - v[j][1];
vr3 = v[i][2] - v[j][2];
vr1 *= dt;
vr2 *= dt;
vr3 *= dt;
// normal component
vnnr = vr1*delx + vr2*dely + vr3*delz;
vn1 = delx*vnnr / rsq;
vn2 = dely*vnnr / rsq;
vn3 = delz*vnnr / rsq;
vn1 = delx*vnnr * rsqinv;
vn2 = dely*vnnr * rsqinv;
vn3 = delz*vnnr * rsqinv;
// tangential component
@ -117,22 +115,17 @@ void PairGranNoHistory::compute(int eflag, int vflag)
// relative rotational velocity
wr1 = radi*omega[i][0] + radj*omega[j][0];
wr2 = radi*omega[i][1] + radj*omega[j][1];
wr3 = radi*omega[i][2] + radj*omega[j][2];
wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
wr1 *= dt/r;
wr2 *= dt/r;
wr3 *= dt/r;
// normal damping term
// this definition of DAMP includes the extra 1/r term
// normal forces = Hookian contact + normal velocity damping
xmeff = rmass[i]*rmass[j] / (rmass[i]+rmass[j]);
if (mask[i] & freeze_group_bit) xmeff = rmass[j];
if (mask[j] & freeze_group_bit) xmeff = rmass[i];
damp = xmeff*gamman_dl*vnnr/rsq;
ccel = xkk*(radsum-r)/r - damp;
damp = xmeff*gamman*vnnr*rsqinv;
ccel = xkk*(radsum-r)*rinv - damp;
// relative velocities
@ -145,11 +138,11 @@ void PairGranNoHistory::compute(int eflag, int vflag)
// force normalization
fn = xmu * fabs(ccel*r);
fs = xmeff*gammas_dl*vrel;
fs = xmeff*gammas*vrel;
if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
else ft = 0.0;
// shear friction forces
// tangential force due to tangential velocity damping
fs1 = -ft*vtr1;
fs2 = -ft*vtr2;
@ -164,7 +157,6 @@ void PairGranNoHistory::compute(int eflag, int vflag)
f[i][1] += fy;
f[i][2] += fz;
rinv = 1/r;
tor1 = rinv * (dely*fs3 - delz*fs2);
tor2 = rinv * (delz*fs1 - delx*fs3);
tor3 = rinv * (delx*fs2 - dely*fs1);

8
src/MOLECULE/Install.csh
View File

@ -34,6 +34,8 @@ if ($1 == 1) then
cp dihedral_multi_harmonic.cpp ..
cp dihedral_opls.cpp ..
cp dump_bond.cpp ..
cp fix_bond_break.cpp ..
cp fix_bond_create.cpp ..
cp fix_bond_swap.cpp ..
cp improper.cpp ..
cp improper_cvff.cpp ..
@ -69,6 +71,8 @@ if ($1 == 1) then
cp dihedral_multi_harmonic.h ..
cp dihedral_opls.h ..
cp dump_bond.h ..
cp fix_bond_break.h ..
cp fix_bond_create.h ..
cp fix_bond_swap.h ..
# cp improper.h ..
cp improper_cvff.h ..
@ -109,6 +113,8 @@ else if ($1 == 0) then
rm ../dihedral_multi_harmonic.cpp
rm ../dihedral_opls.cpp
rm ../dump_bond.cpp
rm ../fix_bond_break.cpp
rm ../fix_bond_create.cpp
rm ../fix_bond_swap.cpp
rm ../improper.cpp
rm ../improper_cvff.cpp
@ -144,6 +150,8 @@ else if ($1 == 0) then
rm ../dihedral_multi_harmonic.h
rm ../dihedral_opls.h
rm ../dump_bond.h
rm ../fix_bond_break.h
rm ../fix_bond_create.h
rm ../fix_bond_swap.h
# rm ../improper.h
rm ../improper_cvff.h

30
src/MOLECULE/angle_hybrid.cpp
View File

@ -103,7 +103,6 @@ void AngleHybrid::compute(int eflag, int vflag)
else evflag = 0;
for (m = 0; m < nstyles; m++) {
if (styles[m] == NULL) continue;
neighbor->nanglelist = nanglelist[m];
neighbor->anglelist = anglelist[m];
@ -168,6 +167,8 @@ void AngleHybrid::settings(int narg, char **arg)
error->all("Angle style hybrid cannot use same angle style twice");
if (strcmp(arg[m],"hybrid") == 0)
error->all("Angle style hybrid cannot have hybrid as an argument");
if (strcmp(arg[m],"none") == 0)
error->all("Angle style hybrid cannot have none as an argument");
styles[m] = force->new_angle(arg[m]);
keywords[m] = new char[strlen(arg[m])+1];
strcpy(keywords[m],arg[m]);
@ -185,12 +186,18 @@ void AngleHybrid::coeff(int which, int narg, char **arg)
int ilo,ihi;
force->bounds(arg[0],atom->nangletypes,ilo,ihi);
// 2nd arg = angle style name (harmonic, etc)
// 2nd arg = angle sub-style name
// allow for "none" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
if (m == nstyles) error->all("Angle coeff for hybrid has invalid style");
int none = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else error->all("Angle coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
@ -199,15 +206,19 @@ void AngleHybrid::coeff(int which, int narg, char **arg)
// invoke sub-style coeff() starting with 1st arg
if (styles[m]) styles[m]->coeff(which,narg-1,&arg[1]);
if (!none) styles[m]->coeff(which,narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is NULL for "none", still set setflag
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
if (none) {
setflag[i] = 1;
map[i] = -1;
} else {
setflag[i] = styles[m]->setflag[i];
map[i] = m;
if (styles[m] == NULL) setflag[i] = 1;
else setflag[i] = styles[m]->setflag[i];
}
}
}
@ -217,6 +228,7 @@ void AngleHybrid::coeff(int which, int narg, char **arg)
double AngleHybrid::equilibrium_angle(int i)
{
if (map[i] < 0) error->one("Invoked angle equil angle on angle style none");
return styles[map[i]]->equilibrium_angle(i);
}
@ -265,8 +277,8 @@ void AngleHybrid::read_restart(FILE *fp)
double AngleHybrid::single(int type, int i1, int i2, int i3)
{
if (styles[map[type]]) return styles[map[type]]->single(type,i1,i2,i3);
return 0.0;
if (map[type] < 0) error->one("Invoked angle single on angle style none");
return styles[map[type]]->single(type,i1,i2,i3);
}
/* ----------------------------------------------------------------------

25
src/MOLECULE/bond_hybrid.cpp
View File

@ -102,7 +102,6 @@ void BondHybrid::compute(int eflag, int vflag)
else evflag = 0;
for (m = 0; m < nstyles; m++) {
if (styles[m] == NULL) continue;
neighbor->nbondlist = nbondlist[m];
neighbor->bondlist = bondlist[m];
@ -167,6 +166,8 @@ void BondHybrid::settings(int narg, char **arg)
error->all("Bond style hybrid cannot use same bond style twice");
if (strcmp(arg[m],"hybrid") == 0)
error->all("Bond style hybrid cannot have hybrid as an argument");
if (strcmp(arg[m],"none") == 0)
error->all("Bond style hybrid cannot have none as an argument");
styles[m] = force->new_bond(arg[m]);
keywords[m] = new char[strlen(arg[m])+1];
strcpy(keywords[m],arg[m]);
@ -184,12 +185,18 @@ void BondHybrid::coeff(int narg, char **arg)
int ilo,ihi;
force->bounds(arg[0],atom->nbondtypes,ilo,ihi);
// 2nd arg = bond style name (harmonic, fene, etc)
// 2nd arg = bond sub-style name
// allow for "none" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
if (m == nstyles) error->all("Bond coeff for hybrid has invalid style");
int none = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else error->all("Bond coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
@ -198,14 +205,15 @@ void BondHybrid::coeff(int narg, char **arg)
// invoke sub-style coeff() starting with 1st arg
if (styles[m]) styles[m]->coeff(narg-1,&arg[1]);
if (!none) styles[m]->coeff(narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is NULL for "none", still set setflag
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
map[i] = m;
setflag[i] = 1;
if (none) map[i] = -1;
else map[i] = m;
}
}
@ -223,6 +231,7 @@ void BondHybrid::init_style()
double BondHybrid::equilibrium_distance(int i)
{
if (map[i] < 0) error->one("Invoked bond equil distance on bond style none");
return styles[map[i]]->equilibrium_distance(i);
}
@ -271,8 +280,8 @@ void BondHybrid::read_restart(FILE *fp)
double BondHybrid::single(int type, double rsq, int i, int j)
{
if (styles[map[type]]) return styles[map[type]]->single(type,rsq,i,j);
return 0.0;
if (map[type] < 0) error->one("Invoked bond single on bond style none");
return styles[map[type]]->single(type,rsq,i,j);
}
/* ----------------------------------------------------------------------

2
src/MOLECULE/bond_quartic.cpp
View File

@ -245,7 +245,7 @@ void BondQuartic::init_style()
if (force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 ||
force->special_lj[3] != 1.0)
error->all("Must use special bonds = 1,1,1 with bond style quartic");
error->all("Bond style quartic requires special_bonds = 1,1,1");
}
/* ----------------------------------------------------------------------

25
src/MOLECULE/dihedral_hybrid.cpp
View File

@ -105,7 +105,6 @@ void DihedralHybrid::compute(int eflag, int vflag)
else evflag = 0;
for (m = 0; m < nstyles; m++) {
if (styles[m] == NULL) continue;
neighbor->ndihedrallist = ndihedrallist[m];
neighbor->dihedrallist = dihedrallist[m];
@ -170,6 +169,8 @@ void DihedralHybrid::settings(int narg, char **arg)
error->all("Dihedral style hybrid cannot use same dihedral style twice");
if (strcmp(arg[m],"hybrid") == 0)
error->all("Dihedral style hybrid cannot have hybrid as an argument");
if (strcmp(arg[m],"none") == 0)
error->all("Dihedral style hybrid cannot have none as an argument");
styles[m] = force->new_dihedral(arg[m]);
keywords[m] = new char[strlen(arg[m])+1];
strcpy(keywords[m],arg[m]);
@ -187,12 +188,18 @@ void DihedralHybrid::coeff(int which, int narg, char **arg)
int ilo,ihi;
force->bounds(arg[0],atom->ndihedraltypes,ilo,ihi);
// 2nd arg = dihedral style name (harmonic, etc)
// 2nd arg = dihedral sub-style name
// allow for "none" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
if (m == nstyles) error->all("Dihedral coeff for hybrid has invalid style");
int none = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else error->all("Dihedral coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
@ -201,15 +208,19 @@ void DihedralHybrid::coeff(int which, int narg, char **arg)
// invoke sub-style coeff() starting with 1st arg
if (styles[m]) styles[m]->coeff(which,narg-1,&arg[1]);
if (!none) styles[m]->coeff(which,narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is NULL for "none", still set setflag
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
if (none) {
setflag[i] = 1;
map[i] = -1;
} else {
setflag[i] = styles[m]->setflag[i];
map[i] = m;
if (styles[m] == NULL) setflag[i] = 1;
else setflag[i] = styles[m]->setflag[i];
}
}
}

2
src/MOLECULE/fix_bond_swap.cpp
View File

@ -77,6 +77,8 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
nmax = 0;
alist = NULL;
naccept = foursome = 0;
}
/* ---------------------------------------------------------------------- */

25
src/MOLECULE/improper_hybrid.cpp
View File

@ -105,7 +105,6 @@ void ImproperHybrid::compute(int eflag, int vflag)
else evflag = 0;
for (m = 0; m < nstyles; m++) {
if (styles[m] == NULL) continue;
neighbor->nimproperlist = nimproperlist[m];
neighbor->improperlist = improperlist[m];
@ -170,6 +169,8 @@ void ImproperHybrid::settings(int narg, char **arg)
error->all("Improper style hybrid cannot use same improper style twice");
if (strcmp(arg[m],"hybrid") == 0)
error->all("Improper style hybrid cannot have hybrid as an argument");
if (strcmp(arg[m],"none") == 0)
error->all("Improper style hybrid cannot have none as an argument");
styles[m] = force->new_improper(arg[m]);
keywords[m] = new char[strlen(arg[m])+1];
strcpy(keywords[m],arg[m]);
@ -187,12 +188,18 @@ void ImproperHybrid::coeff(int which, int narg, char **arg)
int ilo,ihi;
force->bounds(arg[0],atom->nimpropertypes,ilo,ihi);
// 2nd arg = improper style name (harmonic, etc)
// 2nd arg = improper sub-style name
// allow for "none" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
if (m == nstyles) error->all("Improper coeff for hybrid has invalid style");
int none = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else error->all("Improper coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
@ -201,15 +208,19 @@ void ImproperHybrid::coeff(int which, int narg, char **arg)
// invoke sub-style coeff() starting with 1st arg
if (styles[m]) styles[m]->coeff(which,narg-1,&arg[1]);
if (!none) styles[m]->coeff(which,narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is NULL for "none", still set setflag
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
if (none) {
setflag[i] = 1;
map[i] = -1;
} else {
setflag[i] = styles[m]->setflag[i];
map[i] = m;
if (styles[m] == NULL) setflag[i] = 1;
else setflag[i] = styles[m]->setflag[i];
}
}
}

4
src/MOLECULE/style_molecule.h
View File

@ -90,10 +90,14 @@ DumpStyle(bond,DumpBond)
#endif
#ifdef FixInclude
#include "fix_bond_break.h"
#include "fix_bond_create.h"
#include "fix_bond_swap.h"
#endif
#ifdef FixClass
FixStyle(bond/break,FixBondBreak)
FixStyle(bond/create,FixBondCreate)
FixStyle(bond/swap,FixBondSwap)
#endif

3
src/PERI/compute_damage_atom.cpp
View File

@ -18,6 +18,7 @@
#include "string.h"
#include "compute_damage_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
@ -72,6 +73,8 @@ void ComputeDamageAtom::init()
void ComputeDamageAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow damage array if necessary
if (atom->nlocal > nmax) {

5
src/USER-ACKLAND/compute_ackland_atom.cpp
View File

@ -19,8 +19,8 @@
#include "string.h"
#include "compute_ackland_atom.h"
#include "atom.h"
#include "modify.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
@ -34,7 +34,6 @@
using namespace LAMMPS_NS;
enum{UNKNOWN,BCC,FCC,HCP,ICO};
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
@ -102,7 +101,7 @@ void ComputeAcklandAtom::compute_peratom()
int *ilist,*jlist,*numneigh,**firstneigh;
int chi[8];
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
// grow structure array if necessary

1
src/atom.cpp
View File

@ -51,6 +51,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
nbonds = nangles = ndihedrals = nimpropers = 0;
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
extra_bond_per_atom = 0;
firstgroupname = NULL;

6
src/atom.h
View File

@ -34,6 +34,7 @@ class Atom : protected Pointers {
int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
int nbonds,nangles,ndihedrals,nimpropers;
int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
int extra_bond_per_atom;
int firstgroup; // store atoms in this group first, -1 if unset
int nfirst; // # of atoms in first group on this proc
@ -51,8 +52,9 @@ class Atom : protected Pointers {
double *radius,*density,*rmass,*vfrac,*s0;
double **x0;
int maxspecial;
int **nspecial,**special;
int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs
int **special; // IDs of 1-2,1-3,1-4 neighs of each atom
int maxspecial; // special[nlocal][maxspecial]
int *num_bond;
int **bond_type;

12
src/compute.cpp
View File

@ -66,9 +66,10 @@ Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
id_pre = NULL;
timeflag = 0;
invoked = 0;
comm_forward = comm_reverse = 0;
invoked_scalar = invoked_vector = invoked_peratom = -1;
// set modify defaults
extra_dof = domain->dimension;
@ -175,3 +176,12 @@ int Compute::matchstep(int ntimestep)
}
return 0;
}
/* ----------------------------------------------------------------------
clean out list of timesteps to call the compute on
------------------------------------------------------------------------- */
void Compute::clearstep()
{
ntime = 0;
}

6
src/compute.h
View File

@ -54,7 +54,10 @@ class Compute : protected Pointers {
int maxtime; // max # of entries time list can hold
int *tlist; // time list of steps the Compute is called on
int invoked; // set when Compute is invoked, to avoid re-invoking
int invoked_flag; // 1 if invoked or accessed this step, 0 if not
int invoked_scalar; // last timestep on which compute_scalar() was invoked
int invoked_vector; // ditto for compute_vector()
int invoked_peratom; // ditto for compute_peratom()
double dof; // degrees-of-freedom for temperature
@ -85,6 +88,7 @@ class Compute : protected Pointers {
void addstep(int);
int matchstep(int);
void clearstep();
virtual double memory_usage() {return 0.0;}

9
src/compute_centro_atom.cpp
View File

@ -14,8 +14,8 @@
#include "string.h"
#include "compute_centro_atom.h"
#include "atom.h"
#include "modify.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
@ -27,8 +27,6 @@
using namespace LAMMPS_NS;
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
ComputeCentroAtom::ComputeCentroAtom(LAMMPS *lmp, int narg, char **arg) :
@ -59,6 +57,9 @@ ComputeCentroAtom::~ComputeCentroAtom()
void ComputeCentroAtom::init()
{
if (force->pair == NULL)
error->all("Compute centro/atom requires a pair style be defined");
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"centro/atom") == 0) count++;
@ -91,7 +92,7 @@ void ComputeCentroAtom::compute_peratom()
int *ilist,*jlist,*numneigh,**firstneigh;
double pairs[66];
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
// grow centro array if necessary

5
src/compute_coord_atom.cpp
View File

@ -15,6 +15,7 @@
#include "stdlib.h"
#include "compute_coord_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
@ -27,8 +28,6 @@
using namespace LAMMPS_NS;
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) :
@ -90,7 +89,7 @@ void ComputeCoordAtom::compute_peratom()
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
// grow coordination array if necessary

5
src/compute_displace_atom.cpp
View File

@ -15,6 +15,7 @@
#include "string.h"
#include "compute_displace_atom.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
@ -23,8 +24,6 @@
using namespace LAMMPS_NS;
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
ComputeDisplaceAtom::ComputeDisplaceAtom(LAMMPS *lmp, int narg, char **arg) :
@ -75,7 +74,7 @@ void ComputeDisplaceAtom::init()
void ComputeDisplaceAtom::compute_peratom()
{
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
// grow local displacement array if necessary

5
src/compute_erotate_sphere.cpp
View File

@ -14,6 +14,7 @@
#include "mpi.h"
#include "compute_erotate_sphere.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "group.h"
@ -21,8 +22,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INERTIA 0.4 // moment of inertia for sphere
/* ---------------------------------------------------------------------- */
@ -75,7 +74,7 @@ void ComputeERotateSphere::init()
double ComputeERotateSphere::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **omega = atom->omega;
double *radius = atom->radius;

12
src/compute_group_group.cpp
View File

@ -18,6 +18,7 @@
#include "mpi.h"
#include "compute_group_group.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "neighbor.h"
@ -28,9 +29,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeGroupGroup::ComputeGroupGroup(LAMMPS *lmp, int narg, char **arg) :
@ -62,7 +60,7 @@ ComputeGroupGroup::~ComputeGroupGroup()
void ComputeGroupGroup::init()
{
if (force->pair == NULL)
error->all("Compute group/group requires pair style be defined");
error->all("Compute group/group requires a pair style be defined");
pair = force->pair;
cutsq = force->pair->cutsq;
@ -85,8 +83,7 @@ void ComputeGroupGroup::init_list(int id, NeighList *ptr)
double ComputeGroupGroup::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked |= INVOKED_VECTOR;
invoked_scalar = invoked_vector = update->ntimestep;
interact();
return scalar;
@ -96,8 +93,7 @@ double ComputeGroupGroup::compute_scalar()
void ComputeGroupGroup::compute_vector()
{
invoked |= INVOKED_SCALAR;
invoked |= INVOKED_VECTOR;
invoked_scalar = invoked_vector = update->ntimestep;
interact();
}

5
src/compute_ke.cpp
View File

@ -14,6 +14,7 @@
#include "mpi.h"
#include "compute_ke.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "group.h"
@ -21,8 +22,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
/* ---------------------------------------------------------------------- */
ComputeKE::ComputeKE(LAMMPS *lmp, int narg, char **arg) :
@ -45,7 +44,7 @@ void ComputeKE::init()
double ComputeKE::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **v = atom->v;
double *rmass = atom->rmass;

5
src/compute_ke_atom.cpp
View File

@ -14,6 +14,7 @@
#include "string.h"
#include "compute_ke_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "force.h"
@ -22,8 +23,6 @@
using namespace LAMMPS_NS;
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
ComputeKEAtom::ComputeKEAtom(LAMMPS *lmp, int narg, char **arg) :
@ -60,7 +59,7 @@ void ComputeKEAtom::init()
void ComputeKEAtom::compute_peratom()
{
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
// grow ke array if necessary

7
src/compute_pe.cpp
View File

@ -15,6 +15,7 @@
#include "string.h"
#include "compute_pe.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
@ -28,8 +29,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
/* ---------------------------------------------------------------------- */
ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) :
@ -73,7 +72,9 @@ ComputePE::ComputePE(LAMMPS *lmp, int narg, char **arg) :
double ComputePE::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
if (update->eflag_global != invoked_scalar)
error->all("Energy was not tallied on needed timestep");
double one = 0.0;
if (pairflag && force->pair)

7
src/compute_pe_atom.cpp
View File

@ -14,6 +14,7 @@
#include "string.h"
#include "compute_pe_atom.h"
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
@ -26,8 +27,6 @@
using namespace LAMMPS_NS;
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
@ -76,7 +75,9 @@ void ComputePEAtom::compute_peratom()
{
int i;
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
if (update->eflag_atom != invoked_peratom)
error->all("Per-atom energy was not tallied on needed timestep");
// grow local energy array if necessary

17
src/compute_pressure.cpp
View File

@ -16,6 +16,7 @@
#include "stdlib.h"
#include "compute_pressure.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
@ -30,9 +31,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) :
@ -159,13 +157,16 @@ void ComputePressure::init()
double ComputePressure::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
if (update->vflag_global != invoked_scalar)
error->all("Virial was not tallied on needed timestep");
// invoke temperature it it hasn't been already
double t;
if (keflag) {
if (temperature->invoked & INVOKED_SCALAR) t = temperature->scalar;
if (temperature->invoked_scalar == update->ntimestep)
t = temperature->scalar;
else t = temperature->compute_scalar();
}
@ -197,13 +198,15 @@ double ComputePressure::compute_scalar()
void ComputePressure::compute_vector()
{
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
if (update->vflag_global != invoked_vector)
error->all("Virial was not tallied on needed timestep");
// invoke temperature it it hasn't been already
double *ke_tensor;
if (keflag) {
if (!(temperature->invoked & INVOKED_VECTOR))
if (temperature->invoked_vector != update->ntimestep)
temperature->compute_vector();
ke_tensor = temperature->vector;
}

10
src/compute_reduce.cpp
View File

@ -27,10 +27,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INVOKED_PERATOM 4
enum{SUM,MINN,MAXX};
enum{X,V,F,COMPUTE,FIX,VARIABLE};
@ -224,7 +220,7 @@ void ComputeReduce::init()
double ComputeReduce::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double one = compute_one(0);
if (mode == SUM)
@ -240,7 +236,7 @@ double ComputeReduce::compute_scalar()
void ComputeReduce::compute_vector()
{
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
for (int m = 0; m < nvalues; m++) onevec[m] = compute_one(m);
if (mode == SUM)
@ -288,7 +284,7 @@ double ComputeReduce::compute_one(int m)
// invoke compute if not previously invoked
} else if (which[m] == COMPUTE) {
if (!(modify->compute[n]->invoked & INVOKED_PERATOM))
if (modify->compute[n]->invoked_peratom != update->ntimestep)
modify->compute[n]->compute_peratom();
if (j == 0) {

7
src/compute_stress_atom.cpp
View File

@ -15,6 +15,7 @@
#include "string.h"
#include "compute_stress_atom.h"
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
@ -29,8 +30,6 @@
using namespace LAMMPS_NS;
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
ComputeStressAtom::ComputeStressAtom(LAMMPS *lmp, int narg, char **arg) :
@ -85,7 +84,9 @@ void ComputeStressAtom::compute_peratom()
{
int i,j;
invoked |= INVOKED_PERATOM;
invoked_peratom = update->ntimestep;
if (update->vflag_atom != invoked_peratom)
error->all("Per-atom virial was not tallied on needed timestep");
// grow local stress array if necessary

8
src/compute_temp.cpp
View File

@ -15,6 +15,7 @@
#include "string.h"
#include "compute_temp.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
@ -24,9 +25,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeTemp::ComputeTemp(LAMMPS *lmp, int narg, char **arg) :
@ -75,7 +73,7 @@ void ComputeTemp::dof_compute()
double ComputeTemp::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **v = atom->v;
double *mass = atom->mass;
@ -109,7 +107,7 @@ void ComputeTemp::compute_vector()
{
int i;
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
double **v = atom->v;
double *mass = atom->mass;

8
src/compute_temp_com.cpp
View File

@ -16,6 +16,7 @@
#include "string.h"
#include "compute_temp_com.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "group.h"
#include "modify.h"
@ -29,9 +30,6 @@ using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeTempCOM::ComputeTempCOM(LAMMPS *lmp, int narg, char **arg) :
@ -86,7 +84,7 @@ double ComputeTempCOM::compute_scalar()
{
double vthermal[3];
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
if (dynamic) masstotal = group->mass(igroup);
group->vcm(igroup,masstotal,vbias);
@ -125,7 +123,7 @@ void ComputeTempCOM::compute_vector()
int i;
double vthermal[3];
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
if (dynamic) masstotal = group->mass(igroup);
group->vcm(igroup,masstotal,vbias);

7
src/compute_temp_deform.cpp
View File

@ -20,6 +20,7 @@
#include "compute_temp_deform.h"
#include "domain.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "fix.h"
@ -31,8 +32,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp
/* ---------------------------------------------------------------------- */
@ -103,7 +102,7 @@ double ComputeTempDeform::compute_scalar()
{
double lamda[3],vstream[3],vthermal[3];
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **x = atom->x;
double **v = atom->v;
@ -152,7 +151,7 @@ void ComputeTempDeform::compute_vector()
{
double lamda[3],vstream[3],vthermal[3];
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
double **x = atom->x;
double **v = atom->v;

8
src/compute_temp_partial.cpp
View File

@ -15,6 +15,7 @@
#include "stdlib.h"
#include "compute_temp_partial.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
@ -25,9 +26,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeTempPartial::ComputeTempPartial(LAMMPS *lmp, int narg, char **arg) :
@ -95,7 +93,7 @@ int ComputeTempPartial::dof_remove(int i)
double ComputeTempPartial::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **v = atom->v;
double *mass = atom->mass;
@ -130,7 +128,7 @@ void ComputeTempPartial::compute_vector()
{
int i;
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
double **v = atom->v;
double *mass = atom->mass;

8
src/compute_temp_ramp.cpp
View File

@ -16,6 +16,7 @@
#include "string.h"
#include "compute_temp_ramp.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "group.h"
#include "modify.h"
@ -30,9 +31,6 @@ using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeTempRamp::ComputeTempRamp(LAMMPS *lmp, int narg, char **arg) :
@ -151,7 +149,7 @@ double ComputeTempRamp::compute_scalar()
{
double fraction,vramp,vthermal[3];
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **x = atom->x;
double **v = atom->v;
@ -193,7 +191,7 @@ void ComputeTempRamp::compute_vector()
int i;
double fraction,vramp,vthermal[3];
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
double **x = atom->x;
double **v = atom->v;

8
src/compute_temp_region.cpp
View File

@ -15,6 +15,7 @@
#include "string.h"
#include "compute_temp_region.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "domain.h"
@ -24,9 +25,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
ComputeTempRegion::ComputeTempRegion(LAMMPS *lmp, int narg, char **arg) :
@ -77,7 +75,7 @@ int ComputeTempRegion::dof_remove(int i)
double ComputeTempRegion::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
double **x = atom->x;
double **v = atom->v;
@ -122,7 +120,7 @@ void ComputeTempRegion::compute_vector()
{
int i;
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
double **x = atom->x;
double **v = atom->v;

13
src/compute_temp_sphere.cpp
View File

@ -15,6 +15,7 @@
#include "string.h"
#include "compute_temp_sphere.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
@ -24,9 +25,6 @@
using namespace LAMMPS_NS;
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INERTIA 0.4 // moment of inertia for sphere
/* ---------------------------------------------------------------------- */
@ -136,11 +134,10 @@ void ComputeTempSphere::dof_compute()
double ComputeTempSphere::compute_scalar()
{
invoked |= INVOKED_SCALAR;
invoked_scalar = update->ntimestep;
if (tempbias) {
if (!(tbias->invoked & INVOKED_SCALAR))
double tmp = tbias->compute_scalar();
if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar();
tbias->remove_bias_all();
}
@ -186,10 +183,10 @@ void ComputeTempSphere::compute_vector()
{
int i;
invoked |= INVOKED_VECTOR;
invoked_vector = update->ntimestep;
if (tempbias) {
if (!(tbias->invoked & INVOKED_VECTOR)) tbias->compute_vector();
if (tbias->invoked_vector != update->ntimestep) tbias->compute_vector();
tbias->remove_bias_all();
}

3
src/delete_atoms.cpp
View File

@ -201,7 +201,10 @@ void DeleteAtoms::delete_overlap(int narg, char **arg)
lmp->init();
// error check on cutoff
// if no pair style, neighbor list will be empty
if (force->pair == NULL)
error->all("Delete_atoms requires a pair style be defined");
if (cut > neighbor->cutneighmax)
error->all("Delete_atoms cutoff > neighbor cutoff");

1
src/delete_bonds.cpp
View File

@ -21,7 +21,6 @@
#include "neighbor.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "group.h"
#include "special.h"
#include "error.h"

24
src/dump_custom.cpp
View File

@ -32,7 +32,7 @@ using namespace LAMMPS_NS;
// customize by adding keyword to 1st enum
enum{TAG,MOL,TYPE,MASS,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ,
enum{ID,MOL,TYPE,MASS,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ,
VX,VY,VZ,FX,FY,FZ,
Q,MUX,MUY,MUZ,RADIUS,OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ,
QUATW,QUATI,QUATJ,QUATK,TQX,TQY,TQZ,
@ -40,8 +40,6 @@ enum{TAG,MOL,TYPE,MASS,X,Y,Z,XS,YS,ZS,XU,YU,ZU,IX,IY,IZ,
enum{LT,LE,GT,GE,EQ,NEQ};
enum{INT,DOUBLE};
#define INVOKED_PERATOM 4 // same as in computes
/* ---------------------------------------------------------------------- */
DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
@ -271,10 +269,14 @@ int DumpCustom::count()
// invoke Computes for per-atom dump quantities
// only if not already invoked
if (ncompute)
for (i = 0; i < ncompute; i++)
if (!(compute[i]->invoked & INVOKED_PERATOM))
if (ncompute) {
int ntimestep = update->ntimestep;
for (i = 0; i < ncompute; i++) {
if (compute[i]->invoked_peratom != ntimestep)
compute[i]->compute_peratom();
compute[i]->invoked_flag = 1;
}
}
// evaluate atom-style Variables for per-atom dump quantities
@ -322,7 +324,7 @@ int DumpCustom::count()
// customize by adding to if statement
if (thresh_array[ithresh] == TAG) {
if (thresh_array[ithresh] == ID) {
int *tag = atom->tag;
for (i = 0; i < nlocal; i++) dchoose[i] = tag[i];
ptr = dchoose;
@ -658,8 +660,8 @@ void DumpCustom::parse_fields(int narg, char **arg)
for (int iarg = 5; iarg < narg; iarg++) {
i = iarg-5;
if (strcmp(arg[iarg],"tag") == 0) {
pack_choice[i] = &DumpCustom::pack_tag;
if (strcmp(arg[iarg],"id") == 0) {
pack_choice[i] = &DumpCustom::pack_id;
vtype[i] = INT;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (!atom->molecule_flag)
@ -1051,7 +1053,7 @@ int DumpCustom::modify_param(int narg, char **arg)
// set keyword type of threshhold
// customize by adding to if statement
if (strcmp(arg[1],"tag") == 0) thresh_array[nthresh] = TAG;
if (strcmp(arg[1],"id") == 0) thresh_array[nthresh] = ID;
else if (strcmp(arg[1],"mol") == 0) thresh_array[nthresh] = MOL;
else if (strcmp(arg[1],"type") == 0) thresh_array[nthresh] = TYPE;
else if (strcmp(arg[1],"mass") == 0) thresh_array[nthresh] = MASS;
@ -1314,7 +1316,7 @@ void DumpCustom::pack_variable(int n)
/* ---------------------------------------------------------------------- */
void DumpCustom::pack_tag(int n)
void DumpCustom::pack_id(int n)
{
int *tag = atom->tag;
int nlocal = atom->nlocal;

2
src/dump_custom.h
View File

@ -88,7 +88,7 @@ class DumpCustom : public Dump {
typedef void (DumpCustom::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_tag(int);
void pack_id(int);
void pack_molecule(int);
void pack_type(int);
void pack_mass(int);

27
src/fix.cpp
View File

@ -53,6 +53,7 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
virial_flag = 0;
no_change_box = 0;
time_integrate = 0;
time_depend = 0;
restart_pbc = 0;
scalar_flag = vector_flag = peratom_flag = 0;
@ -65,17 +66,21 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
// mask settings - same as in modify.cpp
INITIAL_INTEGRATE = 1;
PRE_EXCHANGE = 2;
PRE_NEIGHBOR = 4;
POST_FORCE = 8;
FINAL_INTEGRATE = 16;
END_OF_STEP = 32;
THERMO_ENERGY = 64;
INITIAL_INTEGRATE_RESPA = 128;
POST_FORCE_RESPA = 256;
FINAL_INTEGRATE_RESPA = 512;
MIN_POST_FORCE = 1024;
MIN_ENERGY = 2048;
POST_INTEGRATE = 2;
PRE_EXCHANGE = 4;
PRE_NEIGHBOR = 8;
PRE_FORCE = 16;
POST_FORCE = 32;
FINAL_INTEGRATE = 64;
END_OF_STEP = 128;
THERMO_ENERGY = 256;
INITIAL_INTEGRATE_RESPA = 512;
POST_INTEGRATE_RESPA = 1024;
PRE_FORCE_RESPA = 2048;
POST_FORCE_RESPA = 4096;
FINAL_INTEGRATE_RESPA = 8192;
MIN_POST_FORCE = 16384;
MIN_ENERGY = 32768;
}
/* ---------------------------------------------------------------------- */

15
src/fix.h
View File

@ -34,13 +34,14 @@ class Fix : protected Pointers {
int virial_flag; // 1 if Fix contributes to virial, 0 if not
int no_change_box; // 1 if cannot swap ortho <-> triclinic
int time_integrate; // 1 if fix performs time integration, 0 if no
int time_depend; // 1 if fix is timestep dependent, 0 if not
int restart_pbc; // 1 if fix moves atoms (except integrate)
// so that write_restart must remap to PBC
int scalar_flag; // 0/1 if compute_scalar() function exists
int vector_flag; // 0/1 if compute_vector() function exists
int size_vector; // N = size of global vector
int scalar_vector_freq; // frequency compute s/v data is available at
int scalar_vector_freq; // frequency s/v data is available at
int extscalar; // 0/1 if scalar is intensive/extensive
int extvector; // 0/1/-1 if vector is all int/ext/extlist
int *extlist; // list of 0/1 int/ext for each vec component
@ -57,9 +58,11 @@ class Fix : protected Pointers {
double virial[6]; // accumlated virial
double **vatom; // accumulated per-atom virial
int INITIAL_INTEGRATE,PRE_EXCHANGE,PRE_NEIGHBOR; // mask settings
int POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY;
int INITIAL_INTEGRATE_RESPA,POST_FORCE_RESPA,FINAL_INTEGRATE_RESPA;
int INITIAL_INTEGRATE,POST_INTEGRATE; // mask settings
int PRE_EXCHANGE,PRE_NEIGHBOR;
int PRE_FORCE,POST_FORCE,FINAL_INTEGRATE,END_OF_STEP,THERMO_ENERGY;
int INITIAL_INTEGRATE_RESPA,POST_INTEGRATE_RESPA;
int PRE_FORCE_RESPA,POST_FORCE_RESPA,FINAL_INTEGRATE_RESPA;
int MIN_POST_FORCE,MIN_ENERGY;
Fix(class LAMMPS *, int, char **);
@ -73,8 +76,10 @@ class Fix : protected Pointers {
virtual void setup(int) {}
virtual void min_setup(int) {}
virtual void initial_integrate(int) {}
virtual void post_integrate() {}
virtual void pre_exchange() {}
virtual void pre_neighbor() {}
virtual void pre_force(int) {}
virtual void post_force(int) {}
virtual void final_integrate() {}
virtual void end_of_step() {}
@ -91,6 +96,8 @@ class Fix : protected Pointers {
virtual int maxsize_restart() {return 0;}
virtual void initial_integrate_respa(int, int, int) {}
virtual void post_integrate_respa(int, int) {}
virtual void pre_force_respa(int, int, int) {}
virtual void post_force_respa(int, int, int) {}
virtual void final_integrate_respa(int) {}

17
src/fix_ave_atom.cpp
View File

@ -27,8 +27,6 @@ using namespace LAMMPS_NS;
enum{X,V,F,COMPUTE,FIX,VARIABLE};
#define INVOKED_PERATOM 4 // same as in computes
/* ---------------------------------------------------------------------- */
FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
@ -36,6 +34,8 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 7) error->all("Illegal fix ave/atom command");
time_depend = 1;
nevery = atoi(arg[3]);
nrepeat = atoi(arg[4]);
peratom_freq = atoi(arg[5]);
@ -167,6 +167,7 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(0);
// zero the array since dump may access it on timestep 0
// zero the array since a variable may access it before first run
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
@ -185,9 +186,9 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
nvalid -= (nrepeat-1)*nevery;
if (nvalid < update->ntimestep) nvalid += peratom_freq;
// add nvalid to ALL computes that store invocation times
// add nvalid to all computes that store invocation times
// since don't know a priori which are invoked by this fix
// once in end_of_step() can just set timestep for ones actually invoked
// once in end_of_step() can set timestep for ones actually invoked
modify->addstep_compute_all(nvalid);
}
@ -268,7 +269,8 @@ void FixAveAtom::end_of_step()
// skip if not step which requires doing something
if (update->ntimestep != nvalid) return;
int ntimestep = update->ntimestep;
if (ntimestep != nvalid) return;
// zero if first step
@ -307,7 +309,8 @@ void FixAveAtom::end_of_step()
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[n];
if (!(compute->invoked & INVOKED_PERATOM)) compute->compute_peratom();
if (compute->invoked_peratom != ntimestep) compute->compute_peratom();
compute->invoked_flag = 1;
if (j == 0) {
double *compute_scalar = compute->scalar_atom;
@ -351,7 +354,7 @@ void FixAveAtom::end_of_step()
}
irepeat = 0;
nvalid = update->ntimestep+peratom_freq - (nrepeat-1)*nevery;
nvalid = ntimestep+peratom_freq - (nrepeat-1)*nevery;
modify->addstep_compute(nvalid);
// average the final result for the Nfreq timestep

3
src/fix_ave_force.cpp
View File

@ -41,6 +41,9 @@ FixAveForce::FixAveForce(LAMMPS *lmp, int narg, char **arg) :
else yvalue = atof(arg[4]);
if (strcmp(arg[5],"NULL") == 0) zflag = 0;
else zvalue = atof(arg[5]);
foriginal_all[0] = foriginal_all[1] =
foriginal_all[2] = foriginal_all[3] = 0.0;
}
/* ---------------------------------------------------------------------- */

19
src/fix_ave_spatial.cpp
View File

@ -38,7 +38,6 @@ enum{BOX,LATTICE,REDUCED};
enum{ONE,RUNNING,WINDOW};
#define BIG 1000000000
#define INVOKED_PERATOM 4 // same as in computes
/* ---------------------------------------------------------------------- */
@ -50,6 +49,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
MPI_Comm_rank(world,&me);
no_change_box = 1;
time_depend = 1;
nevery = atoi(arg[3]);
nrepeat = atoi(arg[4]);
@ -150,7 +150,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
// optional args
normflag = ALL;
scaleflag = BOX;
scaleflag = LATTICE;
fp = NULL;
ave = ONE;
@ -315,9 +315,9 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
nvalid -= (nrepeat-1)*nevery;
if (nvalid < update->ntimestep) nvalid += nfreq;
// add nvalid to ALL computes that store invocation times
// add nvalid to all computes that store invocation times
// since don't know a priori which are invoked by this fix
// once in end_of_step() can just set timestep for ones actually invoked
// once in end_of_step() can set timestep for ones actually invoked
modify->addstep_compute_all(nvalid);
}
@ -417,7 +417,8 @@ void FixAveSpatial::end_of_step()
// skip if not step which requires doing something
if (update->ntimestep != nvalid) return;
int ntimestep = update->ntimestep;
if (ntimestep != nvalid) return;
// if computing the first sample, setup layers
// compute current origin = boundary for some layer
@ -621,7 +622,8 @@ void FixAveSpatial::end_of_step()
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[n];
if (!(compute->invoked & INVOKED_PERATOM)) compute->compute_peratom();
if (compute->invoked_peratom != ntimestep) compute->compute_peratom();
compute->invoked_flag = 1;
double *scalar = compute->scalar_atom;
double **vector = compute->vector_atom;
int jm1 = j - 1;
@ -699,7 +701,7 @@ void FixAveSpatial::end_of_step()
}
irepeat = 0;
nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery;
nvalid = ntimestep+nfreq - (nrepeat-1)*nevery;
modify->addstep_compute(nvalid);
// time average across samples
@ -782,7 +784,7 @@ void FixAveSpatial::end_of_step()
// output result to file
if (fp && me == 0) {
fprintf(fp,"%d %d\n",update->ntimestep,nlayers);
fprintf(fp,"%d %d\n",ntimestep,nlayers);
for (m = 0; m < nlayers; m++) {
fprintf(fp," %d %g %g",m+1,coord[m],count_total[m]/norm);
for (i = 0; i < nvalues; i++) fprintf(fp," %g",values_total[m][i]/norm);
@ -803,6 +805,7 @@ double FixAveSpatial::compute_vector(int n)
int ivalue = n % nvalues;
int ilayer = n / nvalues;
if (ilayer >= nlayers) return 0.0;
if (values_total == NULL) return 0.0;
return values_total[ilayer][ivalue]/norm;
}

24
src/fix_ave_time.cpp
View File

@ -32,9 +32,6 @@ using namespace LAMMPS_NS;
enum{COMPUTE,FIX,VARIABLE};
enum{ONE,RUNNING,WINDOW};
#define INVOKED_SCALAR 1 // same as in computes
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
@ -42,6 +39,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 7) error->all("Illegal fix ave/time command");
time_depend = 1;
MPI_Comm_rank(world,&me);
nevery = atoi(arg[3]);
@ -248,9 +247,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
nvalid -= (nrepeat-1)*nevery;
if (nvalid < update->ntimestep) nvalid += nfreq;
// add nvalid to ALL computes that store invocation times
// add nvalid to all computes that store invocation times
// since don't know a priori which are invoked by this fix
// once in end_of_step() can just set timestep for ones actually invoked
// once in end_of_step() can set timestep for ones actually invoked
modify->addstep_compute_all(nvalid);
}
@ -330,7 +329,8 @@ void FixAveTime::end_of_step()
// skip if not step which requires doing something
if (update->ntimestep != nvalid) return;
int ntimestep = update->ntimestep;
if (ntimestep != nvalid) return;
// zero if first step
@ -351,11 +351,13 @@ void FixAveTime::end_of_step()
Compute *compute = modify->compute[m];
if (argindex[i] == 0) {
if (compute->invoked & INVOKED_SCALAR) vector[i] += compute->scalar;
if (compute->invoked_scalar == ntimestep) vector[i] += compute->scalar;
else vector[i] += compute->compute_scalar();
compute->invoked_flag = 1;
} else {
if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector();
if (compute->invoked_vector != ntimestep) compute->compute_vector();
vector[i] += compute->vector[argindex[i]-1];
compute->invoked_flag = 1;
}
// access fix fields, guaranteed to be ready
@ -383,7 +385,7 @@ void FixAveTime::end_of_step()
}
irepeat = 0;
nvalid = update->ntimestep+nfreq - (nrepeat-1)*nevery;
nvalid = ntimestep+nfreq - (nrepeat-1)*nevery;
modify->addstep_compute(nvalid);
// average the final result for the Nfreq timestep
@ -395,7 +397,7 @@ void FixAveTime::end_of_step()
// if ave = RUNNING, combine with all previous Nfreq timestep values
// if ave = WINDOW, combine with nwindow most recent Nfreq timestep values
if (update->ntimestep >= startstep) {
if (ntimestep >= startstep) {
if (ave == ONE) {
for (i = 0; i < nvalues; i++) vector_total[i] = vector[i];
norm = 1;
@ -424,7 +426,7 @@ void FixAveTime::end_of_step()
// output result to file
if (fp && me == 0) {
fprintf(fp,"%d",update->ntimestep);
fprintf(fp,"%d",ntimestep);
for (i = 0; i < nvalues; i++) fprintf(fp," %g",vector_total[i]/norm);
fprintf(fp,"\n");
fflush(fp);

1
src/fix_deform.cpp
View File

@ -46,6 +46,7 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
box_change = 1;
no_change_box = 1;
time_depend = 1;
nevery = atoi(arg[3]);
if (nevery <= 0) error->all("Illegal fix deform command");

1
src/fix_deposit.cpp
View File

@ -37,6 +37,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
if (narg < 7) error->all("Illegal fix deposit command");
restart_global = 1;
time_depend = 1;
// required args

1
src/fix_dt_reset.cpp
View File

@ -43,6 +43,7 @@ FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 7) error->all("Illegal fix dt/reset command");
time_depend = 1;
scalar_flag = 1;
vector_flag = 1;
size_vector = 1;

2
src/fix_gravity.cpp
View File

@ -32,6 +32,8 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 5) error->all("Illegal fix gravity command");
time_depend = 1;
magnitude = atof(arg[3]);
if (strcmp(arg[4],"chute") == 0) {

1
src/fix_indent.cpp
View File

@ -38,6 +38,7 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 4) error->all("Illegal fix indent command");
time_depend = 1;
scalar_flag = 1;
vector_flag = 1;
size_vector = 3;

2
src/fix_langevin.cpp
View File

@ -40,6 +40,8 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 7) error->all("Illegal fix langevin command");
time_depend = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
t_period = atof(arg[5]);

2
src/fix_nve_limit.cpp
View File

@ -37,6 +37,8 @@ FixNVELimit::FixNVELimit(LAMMPS *lmp, int narg, char **arg) :
extscalar = 1;
xlimit = atof(arg[3]);
ncount = 0;
}
/* ---------------------------------------------------------------------- */

6
src/fix_nvt.cpp
View File

@ -62,8 +62,6 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
if (t_period <= 0.0) error->all("Fix nvt period must be > 0.0");
t_freq = 1.0 / t_period;
eta = eta_dot = 0.0;
// create a new compute temp style
// id = fix-ID + temp, compute group = fix group
@ -86,6 +84,10 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
// Nose/Hoover temp init
eta = eta_dot = 0.0;
}
/* ---------------------------------------------------------------------- */

25
src/fix_recenter.cpp
View File

@ -28,6 +28,8 @@
using namespace LAMMPS_NS;
enum{BOX,LATTICE,FRACTION};
/* ---------------------------------------------------------------------- */
FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
@ -52,7 +54,7 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
// optional args
group2bit = groupbit;
scaleflag = 1;
scaleflag = LATTICE;
int iarg = 6;
while (iarg < narg) {
@ -62,9 +64,9 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
group2bit = group->bitmask[igroup2];
iarg += 2;
} else if (strcmp(arg[iarg],"units") == 0) {
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else if (strcmp(arg[iarg+1],"fraction") == 0) scaleflag = 2;
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE;
else if (strcmp(arg[iarg+1],"fraction") == 0) scaleflag = FRACTION;
else error->all("Illegal fix recenter command");
iarg += 2;
} else error->all("Illegal fix recenter command");
@ -72,11 +74,11 @@ FixRecenter::FixRecenter(LAMMPS *lmp, int narg, char **arg) :
// scale xcom,ycom,zcom
if (scaleflag == 1 && domain->lattice == NULL)
if (scaleflag == LATTICE && domain->lattice == NULL)
error->all("Use of fix recenter with undefined lattice");
double xscale,yscale,zscale;
if (scaleflag == 1) {
if (scaleflag == LATTICE) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
@ -140,7 +142,7 @@ void FixRecenter::initial_integrate(int vflag)
double xtarget,ytarget,ztarget;
double *bboxlo,*bboxhi;
if (scaleflag == 2) {
if (scaleflag == FRACTION) {
if (domain->triclinic == 0) {
bboxlo = domain->boxlo;
bboxhi = domain->boxhi;
@ -151,15 +153,18 @@ void FixRecenter::initial_integrate(int vflag)
}
if (xinitflag) xtarget = xinit;
else if (scaleflag == 2) xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]);
else if (scaleflag == FRACTION)
xtarget = bboxlo[0] + xcom*(bboxhi[0] - bboxlo[0]);
else xtarget = xcom;
if (yinitflag) ytarget = yinit;
else if (scaleflag == 2) ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]);
else if (scaleflag == FRACTION)
ytarget = bboxlo[1] + ycom*(bboxhi[1] - bboxlo[1]);
else ytarget = ycom;
if (zinitflag) ztarget = zinit;
else if (scaleflag == 2) ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]);
else if (scaleflag == FRACTION)
ztarget = bboxlo[2] + zcom*(bboxhi[2] - bboxlo[2]);
else ztarget = zcom;
// current COM

2
src/fix_shake.cpp
View File

@ -64,7 +64,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
grow_arrays(atom->nmax);
atom->add_callback(0);
// set comm size needed by this Fix
// set comm size needed by this fix
comm_forward = 3;

2
src/fix_spring_self.cpp
View File

@ -72,6 +72,8 @@ FixSpringSelf::FixSpringSelf(LAMMPS *lmp, int narg, char **arg) :
xoriginal[i][2] = x[i][2] + zbox*zprd;
} else xoriginal[i][0] = xoriginal[i][1] = xoriginal[i][2] = 0.0;
}
espring = 0.0;
}
/* ---------------------------------------------------------------------- */

1
src/fix_wall_lj126.cpp
View File

@ -35,6 +35,7 @@ FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 8) error->all("Illegal fix wall/lj126 command");
time_depend = 1;
scalar_flag = 1;
vector_flag = 1;
size_vector = 3;

1
src/fix_wall_lj93.cpp
View File

@ -31,6 +31,7 @@ FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 8) error->all("Illegal fix wall/lj93 command");
time_depend = 1;
scalar_flag = 1;
vector_flag = 1;
size_vector = 3;

22
src/fix_wall_reflect.cpp
View File

@ -16,7 +16,6 @@
#include "atom.h"
#include "modify.h"
#include "domain.h"
#include "comm.h"
#include "error.h"
using namespace LAMMPS_NS;
@ -52,29 +51,14 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
int FixWallReflect::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= POST_INTEGRATE;
mask |= POST_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixWallReflect::init()
{
// warn if any integrate fix comes after this one
int after = 0;
int flag = 0;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(id,modify->fix[i]->id) == 0) after = 1;
else if ((modify->fmask[i] & INITIAL_INTEGRATE) && after) flag = 1;
}
if (flag && comm->me == 0)
error->warning("Fix wall/reflect should come after all other integration fixes");
}
/* ---------------------------------------------------------------------- */
void FixWallReflect::initial_integrate(int vflag)
void FixWallReflect::post_integrate()
{
double xlo = domain->boxlo[0];
double xhi = domain->boxhi[0];

3
src/fix_wall_reflect.h
View File

@ -22,8 +22,7 @@ class FixWallReflect : public Fix {
public:
FixWallReflect(class LAMMPS *, int, char **);
int setmask();
void init();
void initial_integrate(int);
void post_integrate();
private:
int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag;

2
src/fix_wiggle.cpp
View File

@ -30,6 +30,8 @@ FixWiggle::FixWiggle(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 6) error->all("Illegal fix wiggle command");
time_depend = 1;
// parse command-line args
if (strcmp(arg[3],"x") == 0) axis = 0;

101
src/force.cpp
View File

@ -56,6 +56,7 @@ Force::Force(LAMMPS *lmp) : Pointers(lmp)
special_lj[1] = special_lj[2] = special_lj[3] = 0.0;
special_coul[1] = special_coul[2] = special_coul[3] = 0.0;
special_dihedral = 0;
special_extra = 0;
dielectric = 1.0;
@ -353,54 +354,72 @@ void Force::set_special(int narg, char **arg)
{
if (narg == 0) error->all("Illegal special_bonds command");
if (strcmp(arg[0],"charmm") == 0) {
if (narg != 1) error->all("Illegal special_bonds command");
special_lj[1] = 0.0;
special_lj[2] = 0.0;
special_lj[3] = 0.0;
special_coul[1] = 0.0;
special_coul[2] = 0.0;
special_coul[3] = 0.0;
special_dihedral = 0;
return;
}
if (strcmp(arg[0],"amber") == 0) {
if (narg != 1) error->all("Illegal special_bonds command");
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"amber") == 0) {
if (iarg+1 > narg) error->all("Illegal special_bonds command");
special_lj[1] = 0.0;
special_lj[2] = 0.0;
special_lj[3] = 0.5;
special_coul[1] = 0.0;
special_coul[2] = 0.0;
special_coul[3] = 5.0/6.0;
special_dihedral = 0;
return;
}
int iarg;
if (strcmp(arg[0],"dihedral") == 0) {
special_dihedral = 1;
iarg = 1;
} else if (strcmp(arg[0],"explicit") == 0) {
special_dihedral = 0;
iarg = 1;
} else {
special_dihedral = 0;
iarg = 0;
}
if (narg-iarg == 3) {
special_lj[1] = special_coul[1] = atof(arg[iarg+0]);
special_lj[2] = special_coul[2] = atof(arg[iarg+1]);
special_lj[3] = special_coul[3] = atof(arg[iarg+2]);
} else if (narg-iarg == 6) {
special_lj[1] = atof(arg[iarg+0]);
special_lj[2] = atof(arg[iarg+1]);
special_lj[3] = atof(arg[iarg+2]);
special_coul[1] = atof(arg[iarg+3]);
special_coul[2] = atof(arg[iarg+4]);
special_coul[3] = atof(arg[iarg+5]);
iarg += 1;
} else if (strcmp(arg[iarg],"charmm") == 0) {
if (iarg+1 > narg) error->all("Illegal special_bonds command");
special_lj[1] = 0.0;
special_lj[2] = 0.0;
special_lj[3] = 0.0;
special_coul[1] = 0.0;
special_coul[2] = 0.0;
special_coul[3] = 0.0;
iarg += 1;
} else if (strcmp(arg[iarg],"fene") == 0) {
if (iarg+1 > narg) error->all("Illegal special_bonds command");
special_lj[1] = 0.0;
special_lj[2] = 1.0;
special_lj[3] = 1.0;
special_coul[1] = 0.0;
special_coul[2] = 1.0;
special_coul[3] = 1.0;
iarg += 1;
} else if (strcmp(arg[iarg],"lj/coul") == 0) {
if (iarg+4 > narg) error->all("Illegal special_bonds command");
special_lj[1] = special_coul[1] = atof(arg[iarg+1]);
special_lj[2] = special_coul[2] = atof(arg[iarg+2]);
special_lj[3] = special_coul[3] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"lj") == 0) {
if (iarg+4 > narg) error->all("Illegal special_bonds command");
special_lj[1] = atof(arg[iarg+1]);
special_lj[2] = atof(arg[iarg+2]);
special_lj[3] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"coul") == 0) {
if (iarg+4 > narg) error->all("Illegal special_bonds command");
special_coul[1] = atof(arg[iarg+1]);
special_coul[2] = atof(arg[iarg+2]);
special_coul[3] = atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
if (iarg+2 > narg) error->all("Illegal special_bonds command");
if (strcmp(arg[iarg+1],"no") == 0) special_dihedral = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) special_dihedral = 1;
else error->all("Illegal special_bonds command");
iarg += 2;
} else if (strcmp(arg[iarg],"extra") == 0) {
if (iarg+2 > narg) error->all("Illegal special_bonds command");
special_extra = atoi(arg[iarg+1]);
iarg += 2;
} else error->all("Illegal special_bonds command");
}
for (int i = 1; i <= 3; i++)
if (special_lj[i] < 0.0 || special_lj[i] > 1.0 ||
special_coul[i] < 0.0 || special_coul[i] > 1.0)
error->all("Illegal special_bonds command");
if (special_extra < 0) error->all("Illegal special_bonds command");
}
/* ----------------------------------------------------------------------

1
src/force.h
View File

@ -54,6 +54,7 @@ class Force : protected Pointers {
double special_coul[4]; // 1-2, 1-3, 1-4 prefactors for Coulombics
int special_dihedral; // 0 if defined dihedrals are ignored
// 1 if only weight 1,4 atoms if in a dihedral
int special_extra; // extra space for added bonds
Force(class LAMMPS *);
~Force();

6
src/input.cpp
View File

@ -1080,13 +1080,15 @@ void Input::shape()
void Input::special_bonds()
{
// store 1-3,1-4 and dihedral flag values before change
// store 1-3,1-4 and dihedral/extra flag values before change
// change in 1-2 coeffs will not change the special list
double lj2 = force->special_lj[2];
double lj3 = force->special_lj[3];
double coul2 = force->special_coul[2];
double coul3 = force->special_coul[3];
int dihedral = force->special_dihedral;
int extra = force->special_extra;
force->set_special(narg,arg);
@ -1095,7 +1097,7 @@ void Input::special_bonds()
if (domain->box_exist && atom->molecular) {
if (lj2 != force->special_lj[2] || lj3 != force->special_lj[3] ||
coul2 != force->special_coul[2] || coul3 != force->special_coul[3] ||
dihedral != force->special_dihedral) {
dihedral != force->special_dihedral || extra != force->special_extra) {
Special special(lmp);
special.build();
}

5
src/integrate.cpp
View File

@ -13,6 +13,7 @@
#include "stdlib.h"
#include "integrate.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
@ -104,6 +105,8 @@ void Integrate::ev_set(int ntimestep)
if (elist_atom[i]->matchstep(ntimestep)) break;
if (i < nelist_atom) eflag_atom = 2;
if (eflag_global) update->eflag_global = update->ntimestep;
if (eflag_atom) update->eflag_atom = update->ntimestep;
eflag = eflag_global + eflag_atom;
int vflag_global = 0;
@ -116,5 +119,7 @@ void Integrate::ev_set(int ntimestep)
if (vlist_atom[i]->matchstep(ntimestep)) break;
if (i < nvlist_atom) vflag_atom = 4;
if (vflag_global) update->vflag_global = update->ntimestep;
if (vflag_atom) update->vflag_atom = update->ntimestep;
vflag = vflag_global + vflag_atom;
}

5
src/min.cpp
View File

@ -14,6 +14,7 @@
#include "stdlib.h"
#include "string.h"
#include "min.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "error.h"
@ -116,6 +117,8 @@ void Min::ev_set(int ntimestep)
if (elist_atom[i]->matchstep(ntimestep)) break;
if (i < nelist_atom) eflag_atom = 2;
if (eflag_global) update->eflag_global = update->ntimestep;
if (eflag_atom) update->eflag_atom = update->ntimestep;
eflag = eflag_global + eflag_atom;
int vflag_global = 0;
@ -128,5 +131,7 @@ void Min::ev_set(int ntimestep)
if (vlist_atom[i]->matchstep(ntimestep)) break;
if (i < nvlist_atom) vflag_atom = 4;
if (vflag_global) update->vflag_global = update->ntimestep;
if (vflag_atom) update->vflag_atom = update->ntimestep;
vflag = vflag_global + vflag_atom;
}

4
src/min_cg.cpp
View File

@ -170,8 +170,8 @@ void MinCG::run()
// set output->next values to this timestep
// call eng_force to insure vflag is set when forces computed
// output->write does final output for thermo, dump, restart files
// add ntimestep to ALL computes that store invocation times
// since just hardwired call to thermo/dumps and they may not be ready
// add ntimestep to all computes that store invocation times
// since just hardwired call to thermo/dumps and computes may not be ready
if (niter < update->nsteps) {
niter++;

141
src/modify.cpp
View File

@ -37,36 +37,44 @@ using namespace LAMMPS_NS;
// mask settings - same as in fix.cpp
#define INITIAL_INTEGRATE 1
#define PRE_EXCHANGE 2
#define PRE_NEIGHBOR 4
#define POST_FORCE 8
#define FINAL_INTEGRATE 16
#define END_OF_STEP 32
#define THERMO_ENERGY 64
#define INITIAL_INTEGRATE_RESPA 128
#define POST_FORCE_RESPA 256
#define FINAL_INTEGRATE_RESPA 512
#define MIN_POST_FORCE 1024
#define MIN_ENERGY 2048
#define POST_INTEGRATE 2
#define PRE_EXCHANGE 4
#define PRE_NEIGHBOR 8
#define PRE_FORCE 16
#define POST_FORCE 32
#define FINAL_INTEGRATE 64
#define END_OF_STEP 128
#define THERMO_ENERGY 256
#define INITIAL_INTEGRATE_RESPA 512
#define POST_INTEGRATE_RESPA 1024
#define PRE_FORCE_RESPA 2048
#define POST_FORCE_RESPA 4096
#define FINAL_INTEGRATE_RESPA 8192
#define MIN_POST_FORCE 16384
#define MIN_ENERGY 32768
/* ---------------------------------------------------------------------- */
Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
{
nfix = maxfix = 0;
n_initial_integrate = 0;
n_initial_integrate = n_post_integrate = 0;
n_pre_exchange = n_pre_neighbor = 0;
n_post_force = n_final_integrate = n_end_of_step = n_thermo_energy = 0;
n_initial_integrate_respa = n_post_force_respa = n_final_integrate_respa = 0;
n_pre_force = n_post_force = 0;
n_final_integrate = n_end_of_step = n_thermo_energy = 0;
n_initial_integrate_respa = n_post_integrate_respa = 0;
n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
n_min_post_force = n_min_energy = 0;
fix = NULL;
fmask = NULL;
list_initial_integrate = NULL;
list_initial_integrate = list_post_integrate = NULL;
list_pre_exchange = list_pre_neighbor = NULL;
list_post_force = list_final_integrate = list_end_of_step = NULL;
list_pre_force = list_post_force = NULL;
list_final_integrate = list_end_of_step = NULL;
list_thermo_energy = NULL;
list_initial_integrate_respa = list_post_force_respa = NULL;
list_initial_integrate_respa = list_post_integrate_respa = NULL;
list_pre_force_respa = list_post_force_respa = NULL;
list_final_integrate_respa = NULL;
list_min_post_force = list_min_energy = NULL;
@ -101,13 +109,17 @@ Modify::~Modify()
memory->sfree(compute);
delete [] list_initial_integrate;
delete [] list_post_integrate;
delete [] list_pre_exchange;
delete [] list_pre_neighbor;
delete [] list_pre_force;
delete [] list_post_force;
delete [] list_final_integrate;
delete [] list_end_of_step;
delete [] list_thermo_energy;
delete [] list_initial_integrate_respa;
delete [] list_post_integrate_respa;
delete [] list_pre_force_respa;
delete [] list_post_force_respa;
delete [] list_final_integrate_respa;
delete [] list_min_post_force;
@ -134,8 +146,10 @@ void Modify::init()
// create lists of fixes to call at each stage of run
list_init(INITIAL_INTEGRATE,n_initial_integrate,list_initial_integrate);
list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate);
list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor);
list_init(PRE_FORCE,n_pre_force,list_pre_force);
list_init(POST_FORCE,n_post_force,list_post_force);
list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate);
list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step);
@ -143,8 +157,12 @@ void Modify::init()
list_init(INITIAL_INTEGRATE_RESPA,
n_initial_integrate_respa,list_initial_integrate_respa);
list_init(POST_INTEGRATE_RESPA,
n_post_integrate_respa,list_post_integrate_respa);
list_init(POST_FORCE_RESPA,
n_post_force_respa,list_post_force_respa);
list_init(PRE_FORCE_RESPA,
n_pre_force_respa,list_pre_force_respa);
list_init(FINAL_INTEGRATE_RESPA,
n_final_integrate_respa,list_final_integrate_respa);
@ -161,9 +179,21 @@ void Modify::init()
for (i = 0; i < nfix; i++) fix[i]->init();
// init each compute
// add initial step to ALL computes that store invocation times
// reset their invoked flags
// else may not be invoked on initial timestep of this run
// if were invoked on last step of previous run
// would be bad if state of system changed between runs
// add initial timestep to all computes that store invocation times
// since any of them may be invoked by initial thermo
// do not clear stored invocation times within compute,
// b/c some may be holdovers from previous run, like for ave fixes
for (i = 0; i < ncompute; i++) compute[i]->init();
for (i = 0; i < ncompute; i++) {
compute[i]->init();
compute[i]->invoked_scalar = -1;
compute[i]->invoked_vector = -1;
compute[i]->invoked_peratom = -1;
}
modify->addstep_compute_all(update->ntimestep);
// set global flag if any fix has its restart_pbc flag set
@ -222,6 +252,16 @@ void Modify::initial_integrate(int vflag)
fix[list_initial_integrate[i]]->initial_integrate(vflag);
}
/* ----------------------------------------------------------------------
post_integrate call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::post_integrate()
{
for (int i = 0; i < n_post_integrate; i++)
fix[list_post_integrate[i]]->post_integrate();
}
/* ----------------------------------------------------------------------
pre_exchange call only for relevant fixes
------------------------------------------------------------------------- */
@ -243,7 +283,17 @@ void Modify::pre_neighbor()
}
/* ----------------------------------------------------------------------
force adjustment call only for relevant fixes
pre_force call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::pre_force(int vflag)
{
for (int i = 0; i < n_pre_force; i++)
fix[list_pre_force[i]]->pre_force(vflag);
}
/* ----------------------------------------------------------------------
post_force call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::post_force(int vflag)
@ -300,7 +350,27 @@ void Modify::initial_integrate_respa(int vflag, int ilevel, int flag)
}
/* ----------------------------------------------------------------------
rRESPA force adjustment call only for relevant fixes
rRESPA post_integrate call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::post_integrate_respa(int ilevel, int flag)
{
for (int i = 0; i < n_post_integrate_respa; i++)
fix[list_post_integrate_respa[i]]->post_integrate_respa(ilevel,flag);
}
/* ----------------------------------------------------------------------
rRESPA pre_force call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::pre_force_respa(int vflag, int ilevel, int iloop)
{
for (int i = 0; i < n_pre_force_respa; i++)
fix[list_pre_force_respa[i]]->pre_force_respa(vflag,ilevel,iloop);
}
/* ----------------------------------------------------------------------
rRESPA post_force call only for relevant fixes
------------------------------------------------------------------------- */
void Modify::post_force_respa(int vflag, int ilevel, int iloop)
@ -614,39 +684,40 @@ int Modify::find_compute(char *id)
}
/* ----------------------------------------------------------------------
clear invoked flag of all computes
always called before computes are invoked
loop only over computes that store invocation times
clear their invoked flag for this step
called before computes are invoked
------------------------------------------------------------------------- */
void Modify::clearstep_compute()
{
for (int icompute = 0; icompute < ncompute; icompute++)
compute[icompute]->invoked = 0;
for (int icompute = 0; icompute < n_timeflag; icompute++)
compute[list_timeflag[icompute]]->invoked_flag = 0;
}
/* ----------------------------------------------------------------------
loop only over computes that store invocation times
if it was invoked, schedule the next invocation
always called after computes are invoked
if its invoked flag set on this timestep, schedule the next invocation
called after computes are invoked
------------------------------------------------------------------------- */
void Modify::addstep_compute(int ntimestep)
void Modify::addstep_compute(int newstep)
{
for (int icompute = 0; icompute < n_timeflag; icompute++)
if (compute[list_timeflag[icompute]]->invoked)
compute[list_timeflag[icompute]]->addstep(ntimestep);
if (compute[list_timeflag[icompute]]->invoked_flag)
compute[list_timeflag[icompute]]->addstep(newstep);
}
/* ----------------------------------------------------------------------
loop only all computes
loop over all computes
schedule the next invocation for those that store invocation times
called when not sure what computes will be needed on newstep
------------------------------------------------------------------------- */
void Modify::addstep_compute_all(int ntimestep)
void Modify::addstep_compute_all(int newstep)
{
for (int icompute = 0; icompute < ncompute; icompute++)
if (compute[icompute]->timeflag)
compute[icompute]->addstep(ntimestep);
for (int icompute = 0; icompute < n_timeflag; icompute++)
compute[list_timeflag[icompute]]->addstep(newstep);
}
/* ----------------------------------------------------------------------

23
src/modify.h
View File

@ -22,9 +22,11 @@ namespace LAMMPS_NS {
class Modify : protected Pointers {
public:
int nfix,maxfix;
int n_initial_integrate,n_pre_decide,n_pre_exchange,n_pre_neighbor;
int n_post_force,n_final_integrate,n_end_of_step,n_thermo_energy;
int n_initial_integrate_respa,n_post_force_respa,n_final_integrate_respa;
int n_initial_integrate,n_post_integrate,n_pre_exchange,n_pre_neighbor;
int n_pre_force,n_post_force;
int n_final_integrate,n_end_of_step,n_thermo_energy;
int n_initial_integrate_respa,n_post_integrate_respa;
int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
int n_min_post_force,n_min_energy;
int nfix_restart_peratom;
@ -41,15 +43,19 @@ class Modify : protected Pointers {
void init();
void setup(int);
void initial_integrate(int);
void post_integrate();
void pre_decide();
void pre_exchange();
void pre_neighbor();
void pre_force(int);
void post_force(int);
void final_integrate();
void end_of_step();
double thermo_energy();
void initial_integrate_respa(int,int,int);
void post_integrate_respa(int,int);
void pre_force_respa(int,int,int);
void post_force_respa(int,int,int);
void final_integrate_respa(int);
@ -78,12 +84,13 @@ class Modify : protected Pointers {
double memory_usage();
private:
// lists of fixes to apply at different times
int *list_initial_integrate,*list_pre_decide;
// lists of fixes to apply at different stages of timestep
int *list_initial_integrate,*list_post_integrate;
int *list_pre_exchange,*list_pre_neighbor;
int *list_post_force,*list_final_integrate,*list_end_of_step;
int *list_thermo_energy;
int *list_initial_integrate_respa,*list_post_force_respa;
int *list_pre_force,*list_post_force;
int *list_final_integrate,*list_end_of_step,*list_thermo_energy;
int *list_initial_integrate_respa,*list_post_integrate_respa;
int *list_pre_force_respa,*list_post_force_respa;
int *list_final_integrate_respa;
int *list_min_post_force,*list_min_energy;

2
src/pair.h
View File

@ -46,7 +46,7 @@ class Pair : protected Pointers {
class NeighList *list; // standard neighbor list used by most pairs
class NeighList *listhalf; // half list used by some pairs
class NeighList *listfull; // full list used by some pairs
class NeighList *listgranhistory; // granular history list used by some
class NeighList *listgranhistory; // granular history list used by some pairs
class NeighList *listinner; // rRESPA lists used by some pairs
class NeighList *listmiddle;
class NeighList *listouter;

5
src/pair_hybrid.cpp
View File

@ -193,11 +193,14 @@ void PairHybrid::settings(int narg, char **arg)
// count sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word
// one exception is 1st two args of style "lj/coul", which are non-numeric
// need a better way to skip these exceptions
nstyles = 0;
i = 0;
while (i < narg) {
if (strcmp(arg[i],"table") == 0) i++;
if (strcmp(arg[i],"lj/coul") == 0) i += 2;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
nstyles++;
@ -304,7 +307,7 @@ void PairHybrid::coeff(int narg, char **arg)
}
// set setflag and which type pairs map to which sub-style
// if sub-style is none: set hybrid subflag, wipe out map
// if sub-style is none: set hybrid setflag, wipe out map
// else: set hybrid setflag & map only if substyle setflag is set
// previous mappings are wiped out

87
src/read_data.cpp
View File

@ -38,6 +38,8 @@ using namespace LAMMPS_NS;
#define CHUNK 1024
#define DELTA 4
#define NSECTIONS 22 // change when add to header::section_keywords
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
@ -82,14 +84,17 @@ void ReadData::command(int narg, char **arg)
if (screen) fprintf(screen,"Scanning data file ...\n");
open(arg[0]);
header(0);
scan(&atom->bond_per_atom,&atom->angle_per_atom,
&atom->dihedral_per_atom,&atom->improper_per_atom);
scan(atom->bond_per_atom,atom->angle_per_atom,
atom->dihedral_per_atom,atom->improper_per_atom);
fclose(fp);
atom->bond_per_atom += atom->extra_bond_per_atom;
}
MPI_Bcast(&atom->bond_per_atom,1,MPI_INT,0,world);
MPI_Bcast(&atom->angle_per_atom,1,MPI_INT,0,world);
MPI_Bcast(&atom->dihedral_per_atom,1,MPI_INT,0,world);
MPI_Bcast(&atom->improper_per_atom,1,MPI_INT,0,world);
} else
atom->bond_per_atom = atom->angle_per_atom =
atom->dihedral_per_atom = atom->improper_per_atom = 0;
@ -122,9 +127,9 @@ void ReadData::command(int narg, char **arg)
domain->set_local_box();
// read rest of file in free format
// if add a section keyword, add to header::section_keywords and NSECTIONS
int atomflag = 0;
parse_keyword(1,1);
while (strlen(keyword)) {
if (strcmp(keyword,"Atoms") == 0) {
@ -242,7 +247,7 @@ void ReadData::command(int narg, char **arg)
} else {
char str[128];
sprintf(str,"Unknown section in data file: %s",keyword);
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->all(str);
}
@ -281,6 +286,15 @@ void ReadData::header(int flag)
int n;
char *ptr;
char *section_keywords[NSECTIONS] =
{"Atoms","Velocities","Bonds","Angles","Dihedrals","Impropers",
"Masses","Shapes","Dipoles",
"Pair Coeffs","Bond Coeffs","Angle Coeffs",
"Dihedral Coeffs","Improper Coeffs",
"BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs",
"EndBondTorsion Coeffs","AngleTorsion Coeffs",
"AngleAngleTorsion Coeffs","BondBond13 Coeffs","AngleAngle Coeffs"};
// skip 1st line of file
if (me == 0) {
@ -322,6 +336,7 @@ void ReadData::header(int flag)
else if (strstr(line,"angles")) sscanf(line,"%d",&atom->nangles);
else if (strstr(line,"dihedrals")) sscanf(line,"%d",&atom->ndihedrals);
else if (strstr(line,"impropers")) sscanf(line,"%d",&atom->nimpropers);
else if (strstr(line,"atom types")) sscanf(line,"%d",&atom->ntypes);
else if (strstr(line,"bond types")) sscanf(line,"%d",&atom->nbondtypes);
else if (strstr(line,"angle types")) sscanf(line,"%d",&atom->nangletypes);
@ -329,6 +344,10 @@ void ReadData::header(int flag)
sscanf(line,"%d",&atom->ndihedraltypes);
else if (strstr(line,"improper types"))
sscanf(line,"%d",&atom->nimpropertypes);
else if (strstr(line,"extra bond per atom"))
sscanf(line,"%d",&atom->extra_bond_per_atom);
else if (strstr(line,"xlo xhi"))
sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]);
else if (strstr(line,"ylo yhi"))
@ -341,6 +360,17 @@ void ReadData::header(int flag)
} else break;
}
// check that exiting string is a valid section keyword
parse_keyword(1,flag);
for (n = 0; n < NSECTIONS; n++)
if (strcmp(keyword,section_keywords[n]) == 0) break;
if (n == NSECTIONS) {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->all(str);
}
// error check on consistency of header values
if ((atom->nbonds || atom->nbondtypes) &&
@ -914,15 +944,14 @@ void ReadData::impropercoeffs(int which)
proc 0 scans the data file for topology maximums
------------------------------------------------------------------------- */
void ReadData::scan(int *pbond_per, int *pangle_per,
int *pdihedral_per, int *pimproper_per)
void ReadData::scan(int &bond_per_atom, int &angle_per_atom,
int &dihedral_per_atom, int &improper_per_atom)
{
int i,tmp1,tmp2,atom1,atom2,atom3,atom4;
int bond_per,angle_per,dihedral_per,improper_per;
char *eof;
int natoms = static_cast<int> (atom->natoms);
bond_per = angle_per = dihedral_per = improper_per = 0;
bond_per_atom = angle_per_atom = dihedral_per_atom = improper_per_atom = 0;
// allocate topology counting vector
// initially, array length = 1 to natoms
@ -931,8 +960,6 @@ void ReadData::scan(int *pbond_per, int *pangle_per,
int cmax = natoms + 1;
int *count = (int *) memory->smalloc(cmax*sizeof(int),"read_data:count");
parse_keyword(1,0);
while (strlen(keyword)) {
if (strcmp(keyword,"Masses") == 0) skip_lines(atom->ntypes);
@ -1050,9 +1077,9 @@ void ReadData::scan(int *pbond_per, int *pangle_per,
count[atom1]++;
count[atom2]++;
}
for (i = 1; i < cmax; i++) bond_per = MAX(bond_per,count[i]);
if (screen) fprintf(screen," %d = max bonds/atom\n",bond_per);
if (logfile) fprintf(logfile," %d = max bonds/atom\n",bond_per);
for (i = 1; i < cmax; i++) bond_per_atom = MAX(bond_per_atom,count[i]);
if (screen) fprintf(screen," %d = max bonds/atom\n",bond_per_atom);
if (logfile) fprintf(logfile," %d = max bonds/atom\n",bond_per_atom);
} else if (strcmp(keyword,"Angles") == 0) {
@ -1077,9 +1104,9 @@ void ReadData::scan(int *pbond_per, int *pangle_per,
count[atom2]++;
count[atom3]++;
}
for (i = 1; i < cmax; i++) angle_per = MAX(angle_per,count[i]);
if (screen) fprintf(screen," %d = max angles/atom\n",angle_per);
if (logfile) fprintf(logfile," %d = max angles/atom\n",angle_per);
for (i = 1; i < cmax; i++) angle_per_atom = MAX(angle_per_atom,count[i]);
if (screen) fprintf(screen," %d = max angles/atom\n",angle_per_atom);
if (logfile) fprintf(logfile," %d = max angles/atom\n",angle_per_atom);
} else if (strcmp(keyword,"Dihedrals") == 0) {
@ -1108,11 +1135,12 @@ void ReadData::scan(int *pbond_per, int *pangle_per,
count[atom3]++;
count[atom4]++;
}
for (i = 1; i < cmax; i++) dihedral_per = MAX(dihedral_per,count[i]);
for (i = 1; i < cmax; i++)
dihedral_per_atom = MAX(dihedral_per_atom,count[i]);
if (screen)
fprintf(screen," %d = max dihedrals/atom\n",dihedral_per);
fprintf(screen," %d = max dihedrals/atom\n",dihedral_per_atom);
if (logfile)
fprintf(logfile," %d = max dihedrals/atom\n",dihedral_per);
fprintf(logfile," %d = max dihedrals/atom\n",dihedral_per_atom);
} else if (strcmp(keyword,"Impropers") == 0) {
for (i = 1; i < cmax; i++) count[i] = 0;
@ -1140,11 +1168,12 @@ void ReadData::scan(int *pbond_per, int *pangle_per,
count[atom3]++;
count[atom4]++;
}
for (i = 1; i < cmax; i++) improper_per = MAX(improper_per,count[i]);
for (i = 1; i < cmax; i++)
improper_per_atom = MAX(improper_per_atom,count[i]);
if (screen)
fprintf(screen," %d = max impropers/atom\n",improper_per);
fprintf(screen," %d = max impropers/atom\n",improper_per_atom);
if (logfile)
fprintf(logfile," %d = max impropers/atom\n",improper_per);
fprintf(logfile," %d = max impropers/atom\n",improper_per_atom);
} else {
char str[128];
@ -1161,17 +1190,11 @@ void ReadData::scan(int *pbond_per, int *pangle_per,
// error check that topology was specified in file
if ((atom->nbonds && !bond_per) || (atom->nangles && !angle_per) ||
(atom->ndihedrals && !dihedral_per) ||
(atom->nimpropers && !improper_per))
if ((atom->nbonds && !bond_per_atom) ||
(atom->nangles && !angle_per_atom) ||
(atom->ndihedrals && !dihedral_per_atom) ||
(atom->nimpropers && !improper_per_atom))
error->one("Needed topology not in data file");
// return values
*pbond_per = bond_per;
*pangle_per = angle_per;
*pdihedral_per = dihedral_per;
*pimproper_per = improper_per;
}
/* ----------------------------------------------------------------------

2
src/read_data.h
View File

@ -33,7 +33,7 @@ class ReadData : protected Pointers {
char **arg;
void open(char *);
void scan(int *, int *, int *, int*);
void scan(int &, int &, int &, int &);
int reallocate(int **, int, int);
void header(int);
void parse_keyword(int, int);

8
src/respa.cpp
View File

@ -415,6 +415,8 @@ void Respa::recurse(int ilevel)
for (int iloop = 0; iloop < loop[ilevel]; iloop++) {
modify->initial_integrate_respa(vflag,ilevel,0);
if (modify->n_post_integrate_respa)
modify->post_integrate_respa(ilevel,iloop);
if (ilevel) recurse(ilevel-1);
@ -427,6 +429,9 @@ void Respa::recurse(int ilevel)
if (ilevel == nlevels-1) {
modify->initial_integrate_respa(vflag,ilevel,1);
if (modify->n_post_integrate_respa)
modify->post_integrate_respa(ilevel,iloop);
int nflag = neighbor->decide();
if (nflag) {
if (modify->n_pre_exchange) modify->pre_exchange();
@ -446,6 +451,7 @@ void Respa::recurse(int ilevel)
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
}
} else if (ilevel == 0) {
timer->stamp();
comm->communicate();
@ -453,6 +459,8 @@ void Respa::recurse(int ilevel)
}
force_clear(newton[ilevel]);
if (modify->n_pre_force_respa)
modify->pre_force_respa(vflag,ilevel,iloop);
timer->stamp();
if (level_bond == ilevel && force->bond) {

15
src/special.cpp
View File

@ -619,11 +619,13 @@ void Special::combine()
}
// compute global maxspecial, must be at least 1
// add in extra factor from special_bonds command
// allocate correct special array with same nmax, new maxspecial
// previously allocated one must be destroyed
// must force AtomVec class to update its ptr to special
// must make AtomVec class update its ptr to special
MPI_Allreduce(&maxspecial,&atom->maxspecial,1,MPI_INT,MPI_MAX,world);
atom->maxspecial += force->special_extra;
atom->maxspecial = MAX(atom->maxspecial,1);
if (me == 0) {
@ -694,7 +696,7 @@ void Special::combine()
/* ----------------------------------------------------------------------
trim list of 1-4 neighbors by checking defined dihedrals
delete a 1-4 neigh if those atoms are not end points of a defined dihedral
delete a 1-4 neigh if they are not end atoms of a defined dihedral
------------------------------------------------------------------------- */
void Special::dihedral_trim()
@ -827,15 +829,6 @@ void Special::dihedral_trim()
} else for (i = 0; i < nlocal; i++) nspecial[i][2] = nspecial[i][1];
// reset atom->maxspecial
int maxspecial = 0;
for (i = 0; i < nlocal; i++)
maxspecial = MAX(maxspecial,nspecial[i][2]);
MPI_Allreduce(&maxspecial,&atom->maxspecial,1,MPI_INT,MPI_MAX,world);
atom->maxspecial = MAX(atom->maxspecial,1);
// stats on new 1-4 neighbor counts
onefourcount = 0.0;

4
src/style_molecule.h
View File

@ -90,10 +90,14 @@ DumpStyle(bond,DumpBond)
#endif
#ifdef FixInclude
#include "fix_bond_break.h"
#include "fix_bond_create.h"
#include "fix_bond_swap.h"
#endif
#ifdef FixClass
FixStyle(bond/break,FixBondBreak)
FixStyle(bond/create,FixBondCreate)
FixStyle(bond/swap,FixBondSwap)
#endif

371
src/thermo.cpp
View File

@ -60,9 +60,6 @@ enum{INT,FLOAT};
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define INVOKED_SCALAR 1 // same as in computes
#define INVOKED_VECTOR 2
/* ---------------------------------------------------------------------- */
Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
@ -73,8 +70,6 @@ Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
style = new char[n];
strcpy(style,arg[0]);
thermoflag = 1;
// set thermo_modify defaults
normuserflag = 0;
@ -261,6 +256,7 @@ void Thermo::compute(int flag)
int i;
firststep = flag;
int ntimestep = update->ntimestep;
// check for lost atoms
// turn off normflag if natoms = 0 to avoid divide by 0
@ -274,12 +270,13 @@ void Thermo::compute(int flag)
for (i = 0; i < ncompute; i++) {
if (compute_which[i] == 0) {
if (!(computes[i]->invoked & INVOKED_SCALAR))
double tmp = computes[i]->compute_scalar();
if (computes[i]->invoked_scalar != ntimestep)
computes[i]->compute_scalar();
} else {
if (!(computes[i]->invoked & INVOKED_VECTOR))
if (computes[i]->invoked_vector != ntimestep)
computes[i]->compute_vector();
}
computes[i]->invoked_flag = 1;
}
// if lineflag = MULTILINE, prepend step/cpu header line
@ -289,7 +286,7 @@ void Thermo::compute(int flag)
double cpu;
if (flag) cpu = timer->elapsed(TIME_LOOP);
else cpu = 0.0;
loc = sprintf(&line[loc],format_multi,update->ntimestep,cpu);
loc = sprintf(&line[loc],format_multi,ntimestep,cpu);
}
// add each thermo value to line with its specific format
@ -301,7 +298,7 @@ void Thermo::compute(int flag)
}
// kludge for RedStorm timing issue
// if (update->ntimestep == 100) return;
// if (ntimestep == 100) return;
// print line to screen and logfile
@ -819,89 +816,203 @@ void Thermo::create_compute(char *id, char *cstyle, char *extra)
}
/* ----------------------------------------------------------------------
compute a single thermodyanmic value
word is any supported keyword in custom list
compute a single thermodyanmic value, word is any keyword in custom list
called when a variable is evaluated by Variable class
return value as double in answer
return 0 if OK, 1 if str is invalid
customize a new keyword by adding to if statement and error tests
return 0 if str is recoginzed keyword, 1 if unrecognized
customize a new keyword by adding to if statement
------------------------------------------------------------------------- */
int Thermo::evaluate_keyword(char *word, double *answer)
{
// check if Compute pointers exist for keywords that need them
// error if thermo style does not use the Compute
if (strcmp(word,"temp") == 0 || strcmp(word,"press") == 0 ||
strcmp(word,"ke") == 0 || strcmp(word,"etotal") == 0 ||
strcmp(word,"pxx") == 0 || strcmp(word,"pyy") == 0 ||
strcmp(word,"pzz") == 0 || strcmp(word,"pxy") == 0 ||
strcmp(word,"pxz") == 0 || strcmp(word,"pyz") == 0)
if (!temperature)
error->all("Variable uses compute via thermo keyword that thermo does not");
if (strcmp(word,"press") == 0 ||
strcmp(word,"pxx") == 0 || strcmp(word,"pyy") == 0 ||
strcmp(word,"pzz") == 0 || strcmp(word,"pxy") == 0 ||
strcmp(word,"pxz") == 0 || strcmp(word,"pyz") == 0)
if (!pressure)
error->all("Variable uses compute via thermo keyword that thermo does not");
if (strcmp(word,"pe") == 0 || strcmp(word,"etotal") == 0 ||
strcmp(word,"enthalpy") == 0)
if (!pe)
error->all("Variable uses compute via thermo keyword that thermo does not");
// set compute_pe invocation flag for keywords that use energy
// but don't call compute_pe explicitly
if (strcmp(word,"evdwl") == 0 || strcmp(word,"ecoul") == 0 ||
strcmp(word,"epair") == 0 || strcmp(word,"ebond") == 0 ||
strcmp(word,"eangle") == 0 || strcmp(word,"edihed") == 0 ||
strcmp(word,"eimp") == 0 || strcmp(word,"emol") == 0 ||
strcmp(word,"elong") == 0 || strcmp(word,"etail") == 0) {
if (!pe)
error->all("Variable uses compute via thermo keyword that thermo does not");
pe->invoked |= INVOKED_SCALAR;
}
// toggle thermoflag off/on
// so individual compute routines know they are not being called from
// Thermo::compute() which pre-calls Computes
thermoflag = 0;
// invoke the lo-level thermo routine to compute the variable value
// invoke a lo-level thermo routine to compute the variable value
// if keyword requires a compute, is error if thermo doesn't use the compute
// if in middle of run and needed compute is not current, invoke it
// if inbetween runs and needed compute is not current, error
// set invoked flag for pe and pressure for keywords that use them
// this insures tallying on future needed steps via clearstep/addstep
// for keywords that use pe indirectly (evdwl, ebond, etc):
// check if energy was tallied on this timestep and set pe->invoked_flag
if (strcmp(word,"step") == 0) {
compute_step();
dvalue = ivalue;
} else if (strcmp(word,"atoms") == 0) {
compute_atoms();
dvalue = ivalue;
} else if (strcmp(word,"cpu") == 0) {
if (update->whichflag < 0) firststep = 0;
compute_cpu();
} else if (strcmp(word,"temp") == 0) {
if (!temperature)
error->all("Thermo keyword in variable requires thermo to use/init temp");
if (temperature->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else temperature->compute_scalar();
}
else if (strcmp(word,"cpu") == 0) compute_cpu();
compute_temp();
else if (strcmp(word,"temp") == 0) compute_temp();
} else if (strcmp(word,"press") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_scalar();
}
pressure->invoked_flag = 1;
compute_press();
else if (strcmp(word,"press") == 0) compute_press();
else if (strcmp(word,"pe") == 0) compute_pe();
else if (strcmp(word,"ke") == 0) compute_ke();
else if (strcmp(word,"etotal") == 0) compute_etotal();
else if (strcmp(word,"enthalpy") == 0) compute_enthalpy();
} else if (strcmp(word,"pe") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (pe->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pe->compute_scalar();
}
pe->invoked_flag = 1;
compute_pe();
else if (strcmp(word,"evdwl") == 0) compute_evdwl();
else if (strcmp(word,"ecoul") == 0) compute_ecoul();
else if (strcmp(word,"epair") == 0) compute_epair();
else if (strcmp(word,"ebond") == 0) compute_ebond();
else if (strcmp(word,"eangle") == 0) compute_eangle();
else if (strcmp(word,"edihed") == 0) compute_edihed();
else if (strcmp(word,"eimp") == 0) compute_eimp();
else if (strcmp(word,"emol") == 0) compute_emol();
else if (strcmp(word,"elong") == 0) compute_elong();
else if (strcmp(word,"etail") == 0) compute_etail();
} else if (strcmp(word,"ke") == 0) {
if (!temperature)
error->all("Thermo keyword in variable requires thermo to use/init temp");
if (temperature->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else temperature->compute_scalar();
}
compute_ke();
else if (strcmp(word,"vol") == 0) compute_vol();
} else if (strcmp(word,"etotal") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (pe->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pe->compute_scalar();
}
if (!temperature)
error->all("Thermo keyword in variable requires thermo to use/init temp");
if (temperature->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else temperature->compute_scalar();
}
pe->invoked_flag = 1;
compute_etotal();
} else if (strcmp(word,"enthalpy") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (pe->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pe->compute_scalar();
}
if (!temperature)
error->all("Thermo keyword in variable requires thermo to use/init temp");
if (temperature->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else temperature->compute_scalar();
}
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_scalar();
}
pe->invoked_flag = 1;
pressure->invoked_flag = 1;
compute_enthalpy();
} else if (strcmp(word,"evdwl") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_evdwl();
} else if (strcmp(word,"ecoul") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_ecoul();
} else if (strcmp(word,"epair") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_epair();
} else if (strcmp(word,"ebond") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_ebond();
} else if (strcmp(word,"eangle") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_eangle();
} else if (strcmp(word,"edihed") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_edihed();
} else if (strcmp(word,"eimp") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_eimp();
} else if (strcmp(word,"emol") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_emol();
} else if (strcmp(word,"elong") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_elong();
} else if (strcmp(word,"etail") == 0) {
if (!pe)
error->all("Thermo keyword in variable requires thermo to use/init pe");
if (update->eflag_global != update->ntimestep)
error->all("Energy was not tallied on needed timestep");
pe->invoked_flag = 1;
compute_etail();
} else if (strcmp(word,"vol") == 0) compute_vol();
else if (strcmp(word,"lx") == 0) compute_lx();
else if (strcmp(word,"ly") == 0) compute_ly();
else if (strcmp(word,"lz") == 0) compute_lz();
@ -917,16 +1028,73 @@ int Thermo::evaluate_keyword(char *word, double *answer)
else if (strcmp(word,"xz") == 0) compute_xz();
else if (strcmp(word,"yz") == 0) compute_yz();
else if (strcmp(word,"pxx") == 0) compute_pxx();
else if (strcmp(word,"pyy") == 0) compute_pyy();
else if (strcmp(word,"pzz") == 0) compute_pzz();
else if (strcmp(word,"pxy") == 0) compute_pxy();
else if (strcmp(word,"pxz") == 0) compute_pxz();
else if (strcmp(word,"pyz") == 0) compute_pyz();
else if (strcmp(word,"pxx") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_vector();
}
pressure->invoked_flag = 1;
compute_pxx();
else return 1;
} else if (strcmp(word,"pyy") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_vector();
}
pressure->invoked_flag = 1;
compute_pyy();
thermoflag = 1;
} else if (strcmp(word,"pzz") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_vector();
}
pressure->invoked_flag = 1;
compute_pzz();
} else if (strcmp(word,"pxy") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_vector();
}
pressure->invoked_flag = 1;
compute_pxy();
} else if (strcmp(word,"pxz") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_vector();
}
pressure->invoked_flag = 1;
compute_pxz();
} else if (strcmp(word,"pyz") == 0) {
if (!pressure)
error->all("Thermo keyword in variable requires thermo to use/init press");
if (pressure->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable thermo keyword is not current");
else pressure->compute_vector();
}
pressure->invoked_flag = 1;
compute_pyz();
} else return 1;
*answer = dvalue;
return 0;
@ -934,8 +1102,6 @@ int Thermo::evaluate_keyword(char *word, double *answer)
/* ----------------------------------------------------------------------
extraction of Compute, Fix, Variable results
ignore thermoflag:
3 methods only called by Thermo::compute(), not by variable evaluation
compute/fix are normalized by atoms if returning extensive value
variable value is not normalized (formula should normalize if desired)
------------------------------------------------------------------------- */
@ -986,8 +1152,8 @@ void Thermo::compute_variable()
/* ----------------------------------------------------------------------
one method for every keyword thermo can output
thermoflag = 1 if called by Thermo::compute()
thermoflag = 0 if called by evaluate_keyword() via Variable
called by compute() or evaluate_keyword()
compute will have already been called
set ivalue/dvalue if value is integer/double
customize a new keyword by adding a method
------------------------------------------------------------------------- */
@ -1001,8 +1167,7 @@ void Thermo::compute_step()
void Thermo::compute_atoms()
{
if (thermoflag) ivalue = static_cast<int> (natoms);
else ivalue = static_cast<int> (atom->natoms);
ivalue = static_cast<int> (natoms);
}
/* ---------------------------------------------------------------------- */
@ -1017,27 +1182,21 @@ void Thermo::compute_cpu()
void Thermo::compute_temp()
{
if (thermoflag || temperature->invoked & INVOKED_SCALAR)
dvalue = temperature->scalar;
else dvalue = temperature->compute_scalar();
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_press()
{
if (thermoflag || pressure->invoked & INVOKED_SCALAR)
dvalue = pressure->scalar;
else dvalue = pressure->compute_scalar();
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pe()
{
if (thermoflag || pe->invoked & INVOKED_SCALAR)
dvalue = pe->scalar;
else dvalue = pe->compute_scalar();
if (normflag) dvalue /= natoms;
}
@ -1045,10 +1204,7 @@ void Thermo::compute_pe()
void Thermo::compute_ke()
{
if (thermoflag || temperature->invoked & INVOKED_SCALAR)
dvalue = temperature->scalar;
else dvalue = temperature->compute_scalar();
dvalue *= 0.5 * temperature->dof * force->boltz;
if (normflag) dvalue /= natoms;
}
@ -1058,12 +1214,7 @@ void Thermo::compute_ke()
void Thermo::compute_etotal()
{
compute_pe();
double ke;
if (thermoflag || temperature->invoked & INVOKED_SCALAR)
ke = temperature->scalar;
else ke = temperature->compute_scalar();
double ke = temperature->scalar;
ke *= 0.5 * temperature->dof * force->boltz;
if (normflag) ke /= natoms;
dvalue += ke;
@ -1307,8 +1458,6 @@ void Thermo::compute_yz()
void Thermo::compute_pxx()
{
if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR))
pressure->compute_vector();
dvalue = pressure->vector[0];
}
@ -1316,8 +1465,6 @@ void Thermo::compute_pxx()
void Thermo::compute_pyy()
{
if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR))
pressure->compute_vector();
dvalue = pressure->vector[1];
}
@ -1325,8 +1472,6 @@ void Thermo::compute_pyy()
void Thermo::compute_pzz()
{
if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR))
pressure->compute_vector();
dvalue = pressure->vector[2];
}
@ -1334,8 +1479,6 @@ void Thermo::compute_pzz()
void Thermo::compute_pxy()
{
if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR))
pressure->compute_vector();
dvalue = pressure->vector[3];
}
@ -1343,8 +1486,6 @@ void Thermo::compute_pxy()
void Thermo::compute_pxz()
{
if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR))
pressure->compute_vector();
dvalue = pressure->vector[4];
}
@ -1352,7 +1493,5 @@ void Thermo::compute_pxz()
void Thermo::compute_pyz()
{
if (thermoflag == 0 && !(pressure->invoked & INVOKED_VECTOR))
pressure->compute_vector();
dvalue = pressure->vector[5];
}

1
src/thermo.h
View File

@ -60,7 +60,6 @@ class Thermo : protected Pointers {
int ivalue; // integer value to print
double dvalue,natoms; // dvalue = double value to print
int ifield; // which field in thermo output is being computed
int thermoflag; // 1 when called by compute(), 0 from variable eval
int *field2index; // which compute,fix,variable calcs this field
int *argindex; // index into compute,fix scalar,vector

24
src/update.cpp
View File

@ -16,6 +16,9 @@
#include "update.h"
#include "neighbor.h"
#include "force.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "output.h"
#include "memory.h"
#include "error.h"
@ -42,6 +45,8 @@ Update::Update(LAMMPS *lmp) : Pointers(lmp)
firststep = laststep = 0;
beginstep = endstep = 0;
eflag_global = vflag_global = -1;
unit_style = NULL;
set_units("lj");
@ -207,6 +212,10 @@ void Update::create_minimize(int narg, char **arg)
/* ----------------------------------------------------------------------
reset timestep from input script
do not allow dump files or a restart to be defined
do not allow any timestep-dependent fixes to be defined
reset eflag/vflag global so nothing will think eng/virial are current
reset invoked flags of computes, so nothing will think they are current
clear timestep list of computes that store one
------------------------------------------------------------------------- */
void Update::reset_timestep(int narg, char **arg)
@ -219,6 +228,21 @@ void Update::reset_timestep(int narg, char **arg)
if (output->restart && output->last_restart >= 0)
error->all("Cannot reset timestep with restart file already written");
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->time_depend)
error->all("Cannot reset timestep with a time-dependent fix defined");
eflag_global = vflag_global = -1;
for (int i = 0; i < modify->ncompute; i++) {
modify->compute[i]->invoked_scalar = -1;
modify->compute[i]->invoked_vector = -1;
modify->compute[i]->invoked_peratom = -1;
}
for (int i = 0; i < modify->ncompute; i++)
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
ntimestep = atoi(arg[0]);
if (ntimestep < 0) error->all("Timestep must be >= 0");
}

3
src/update.h
View File

@ -30,6 +30,9 @@ class Update : protected Pointers {
int first_update; // 0 before initial update, 1 after
int max_eval; // max force evaluations for minimizer
int eflag_global,eflag_atom; // timestep global/peratom eng is tallied
int vflag_global,vflag_atom; // ditto for virial
char *unit_style;
class Integrate *integrate;

88
src/variable.cpp
View File

@ -43,10 +43,6 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,
CEIL,FLOOR,ROUND,VALUE,ATOMARRAY,TYPEARRAY};
#define INVOKED_SCALAR 1 // same as in computes
#define INVOKED_VECTOR 2
#define INVOKED_PERATOM 4
/* ---------------------------------------------------------------------- */
Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
@ -643,15 +639,15 @@ double Variable::evaluate(char *str, Tree **tree)
char *id = new char[n];
strcpy(id,&word[2]);
if (update->first_update == 0)
error->all("Compute in variable formula before initial run");
int icompute = modify->find_compute(id);
if (icompute < 0) error->all("Invalid compute ID in variable formula");
Compute *compute = modify->compute[icompute];
compute->invoked_flag = 1;
delete [] id;
if (domain->box_exist == 0)
error->all("Variable evalulation before simulation box is defined");
// parse zero or one or two trailing brackets
// point i beyond last bracket
// nbracket = # of bracket pairs
@ -674,9 +670,12 @@ double Variable::evaluate(char *str, Tree **tree)
if (nbracket == 0 && compute->scalar_flag) {
if (!(compute->invoked & INVOKED_SCALAR))
value1 = compute->compute_scalar();
else value1 = compute->scalar;
if (compute->invoked_scalar != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable is not current");
else compute->compute_scalar();
}
value1 = compute->scalar;
if (tree) {
Tree *newtree = new Tree();
newtree->type = VALUE;
@ -691,7 +690,11 @@ double Variable::evaluate(char *str, Tree **tree)
if (index1 > compute->size_vector)
error->all("Compute vector in variable formula is too small");
if (!(compute->invoked & INVOKED_VECTOR)) compute->compute_vector();
if (compute->invoked_vector != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable is not current");
else compute->compute_vector();
}
value1 = compute->vector[index1-1];
if (tree) {
Tree *newtree = new Tree();
@ -708,8 +711,11 @@ double Variable::evaluate(char *str, Tree **tree)
if (tree == NULL)
error->all("Per-atom compute in equal-style variable formula");
if (!(compute->invoked & INVOKED_PERATOM))
compute->compute_peratom();
if (compute->invoked_peratom != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable is not current");
else compute->compute_peratom();
}
Tree *newtree = new Tree();
newtree->type = ATOMARRAY;
newtree->array = compute->scalar_atom;
@ -722,8 +728,11 @@ double Variable::evaluate(char *str, Tree **tree)
} else if (nbracket == 1 && index1 > 0 &&
compute->peratom_flag && compute->size_peratom == 0) {
if (!(compute->invoked & INVOKED_PERATOM))
compute->compute_peratom();
if (compute->invoked_peratom != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable is not current");
else compute->compute_peratom();
}
peratom2global(1,NULL,compute->scalar_atom,1,index1,
tree,treestack,ntreestack,argstack,nargstack);
@ -736,8 +745,11 @@ double Variable::evaluate(char *str, Tree **tree)
error->all("Per-atom compute in equal-style variable formula");
if (index2 > compute->size_peratom)
error->all("Compute vector in variable formula is too small");
if (!(compute->invoked & INVOKED_PERATOM))
compute->compute_peratom();
if (compute->invoked_peratom != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable is not current");
else compute->compute_peratom();
}
Tree *newtree = new Tree();
newtree->type = ATOMARRAY;
newtree->array = &compute->vector_atom[0][index2-1];
@ -752,8 +764,11 @@ double Variable::evaluate(char *str, Tree **tree)
if (index2 > compute->size_peratom)
error->all("Compute vector in variable formula is too small");
if (!(compute->invoked & INVOKED_PERATOM))
compute->compute_peratom();
if (compute->invoked_peratom != update->ntimestep) {
if (update->whichflag < 0)
error->all("Compute used in variable is not current");
else compute->compute_peratom();
}
peratom2global(1,NULL,&compute->vector_atom[0][index2-1],
compute->size_peratom,index1,
tree,treestack,ntreestack,argstack,nargstack);
@ -769,15 +784,14 @@ double Variable::evaluate(char *str, Tree **tree)
char *id = new char[n];
strcpy(id,&word[2]);
if (update->first_update == 0)
error->all("Fix in variable formula before initial run");
int ifix = modify->find_fix(id);
if (ifix < 0) error->all("Invalid fix ID in variable formula");
Fix *fix = modify->fix[ifix];
delete [] id;
if (domain->box_exist == 0)
error->all("Variable evaluation before simulation box is defined");
// parse zero or one or two trailing brackets
// point i beyond last bracket
// nbracket = # of bracket pairs
@ -800,7 +814,8 @@ double Variable::evaluate(char *str, Tree **tree)
if (nbracket == 0 && fix->scalar_flag) {
if (update->ntimestep % fix->scalar_vector_freq)
if (update->whichflag >= 0 &&
update->ntimestep % fix->scalar_vector_freq)
error->all("Fix in variable not computed at compatible time");
value1 = fix->compute_scalar();
if (tree) {
@ -817,7 +832,8 @@ double Variable::evaluate(char *str, Tree **tree)
if (index1 > fix->size_vector)
error->all("Fix vector in variable formula is too small");
if (update->ntimestep % fix->scalar_vector_freq)
if (update->whichflag >= 0 &&
update->ntimestep % fix->scalar_vector_freq)
error->all("Fix in variable not computed at compatible time");
value1 = fix->compute_vector(index1-1);
if (tree) {
@ -835,7 +851,8 @@ double Variable::evaluate(char *str, Tree **tree)
if (tree == NULL)
error->all("Per-atom fix in equal-style variable formula");
if (update->ntimestep % fix->peratom_freq)
if (update->whichflag >= 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
Tree *newtree = new Tree();
newtree->type = ATOMARRAY;
@ -849,7 +866,8 @@ double Variable::evaluate(char *str, Tree **tree)
} else if (nbracket == 1 && index1 > 0 &&
fix->peratom_flag && fix->size_peratom == 0) {
if (update->ntimestep % fix->peratom_freq)
if (update->whichflag >= 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
peratom2global(1,NULL,fix->scalar_atom,1,index1,
tree,treestack,ntreestack,argstack,nargstack);
@ -863,7 +881,8 @@ double Variable::evaluate(char *str, Tree **tree)
error->all("Per-atom fix in equal-style variable formula");
if (index2 > fix->size_peratom)
error->all("Fix vector in variable formula is too small");
if (update->ntimestep % fix->peratom_freq)
if (update->whichflag >= 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
Tree *newtree = new Tree();
newtree->type = ATOMARRAY;
@ -879,7 +898,8 @@ double Variable::evaluate(char *str, Tree **tree)
if (index2 > fix->size_peratom)
error->all("Fix vector in variable formula is too small");
if (update->ntimestep % fix->peratom_freq)
if (update->whichflag >= 0 &&
update->ntimestep % fix->peratom_freq)
error->all("Fix in variable not computed at compatible time");
peratom2global(1,NULL,&fix->vector_atom[0][index2-1],
fix->size_peratom,index1,
@ -984,6 +1004,9 @@ double Variable::evaluate(char *str, Tree **tree)
i = int_between_brackets(str,i,id,1);
i++;
if (domain->box_exist == 0)
error->all("Variable evaluation before simulation box is defined");
// ID between brackets exists: atom value
// empty brackets: atom vector
@ -998,8 +1021,9 @@ double Variable::evaluate(char *str, Tree **tree)
// ----------------
} else {
if (update->first_update == 0)
error->all("Thermo keyword in variable formula before initial run");
if (domain->box_exist == 0)
error->all("Variable evaluation before simulation box is defined");
int flag = output->thermo->evaluate_keyword(word,&value1);
if (flag) error->all("Invalid thermo keyword in variable formula");
if (tree) {

2
src/verlet.cpp
View File

@ -137,6 +137,7 @@ void Verlet::iterate(int n)
// initial time integration
modify->initial_integrate(vflag);
if (modify->n_post_integrate) modify->post_integrate();
// regular communication vs neighbor list rebuild
@ -168,6 +169,7 @@ void Verlet::iterate(int n)
// force computations
force_clear();
if (modify->n_pre_force) modify->pre_force(vflag);
timer->stamp();