From cea02082079637a6eb6199408b53ef47a23f4aa2 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 22 Jun 2007 16:59:17 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@642 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/atom_vec_ellipsoid.cpp | 8 -- src/ASPHERE/atom_vec_ellipsoid.h | 1 - src/ASPHERE/compute_temp_asphere.cpp | 3 +- src/DIPOLE/compute_temp_dipole.cpp | 3 +- src/GRANULAR/atom_vec_granular.cpp | 64 ++++++------- src/GRANULAR/atom_vec_granular.h | 2 +- src/GRANULAR/fix_nve_gran.cpp | 21 +++-- src/GRANULAR/fix_nve_gran.h | 2 +- src/GRANULAR/fix_pour.cpp | 24 ++--- src/KSPACE/ewald.cpp | 2 +- src/KSPACE/pppm.cpp | 2 +- src/atom.cpp | 6 +- src/atom.h | 4 +- src/comm.cpp | 6 +- src/compute_rotate_dipole.cpp | 4 +- src/compute_rotate_gran.cpp | 3 +- src/compute_temp.cpp | 3 +- src/compute_temp_deform.cpp | 2 +- src/compute_temp_ramp.cpp | 2 +- src/compute_temp_region.cpp | 2 +- src/create_box.cpp | 2 +- src/domain.cpp | 3 +- src/domain.h | 2 +- src/dump_custom.cpp | 24 ++--- src/fix_deform.cpp | 14 ++- src/fix_deposit.cpp | 14 +-- src/fix_enforce2d.cpp | 19 +++- src/fix_nvt.cpp | 2 +- src/force.cpp | 1 - src/force.h | 1 - src/input.cpp | 2 +- src/lattice.cpp | 4 +- src/neighbor.cpp | 30 +++--- src/pair.cpp | 2 +- src/read_data.cpp | 2 +- src/read_restart.cpp | 6 +- src/replicate.cpp | 2 +- src/set.cpp | 133 ++++++++++++++++++--------- src/set.h | 2 +- src/style.h | 2 + src/velocity.cpp | 6 +- src/write_restart.cpp | 2 +- 42 files changed, 245 insertions(+), 194 deletions(-) diff --git a/src/ASPHERE/atom_vec_ellipsoid.cpp b/src/ASPHERE/atom_vec_ellipsoid.cpp index 2f8c868eed..e2c5f40b0a 100755 --- a/src/ASPHERE/atom_vec_ellipsoid.cpp +++ b/src/ASPHERE/atom_vec_ellipsoid.cpp @@ -48,14 +48,6 @@ AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp, int narg, char **arg) : atom->angmom_flag = atom->torque_flag = atom->quat_flag = 1; } -/* ---------------------------------------------------------------------- */ - -void AtomVecEllipsoid::init() -{ - if (force->dimension == 2) - error->all("Cannot use atom style ellipsoid for 2d simulation"); -} - /* ---------------------------------------------------------------------- grow atom arrays n = 0 grows arrays by DELTA diff --git a/src/ASPHERE/atom_vec_ellipsoid.h b/src/ASPHERE/atom_vec_ellipsoid.h index ce0c0bfe27..3b8c0265e5 100755 --- a/src/ASPHERE/atom_vec_ellipsoid.h +++ b/src/ASPHERE/atom_vec_ellipsoid.h @@ -22,7 +22,6 @@ class AtomVecEllipsoid : public AtomVec { public: AtomVecEllipsoid(class LAMMPS *, int, char **); virtual ~AtomVecEllipsoid() {} - void init(); void grow(int); void copy(int, int); virtual int pack_comm(int, int *, double *, int, int *); diff --git a/src/ASPHERE/compute_temp_asphere.cpp b/src/ASPHERE/compute_temp_asphere.cpp index 6604435487..6923ef1204 100755 --- a/src/ASPHERE/compute_temp_asphere.cpp +++ b/src/ASPHERE/compute_temp_asphere.cpp @@ -20,6 +20,7 @@ #include "math_extra.h" #include "atom.h" #include "force.h" +#include "domain.h" #include "modify.h" #include "fix.h" #include "group.h" @@ -73,7 +74,7 @@ void ComputeTempAsphere::init() void ComputeTempAsphere::recount() { double natoms = group->count(igroup); - dof = force->dimension * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; // add rotational degrees of freedom diff --git a/src/DIPOLE/compute_temp_dipole.cpp b/src/DIPOLE/compute_temp_dipole.cpp index b0ed8c267e..dac8f11671 100644 --- a/src/DIPOLE/compute_temp_dipole.cpp +++ b/src/DIPOLE/compute_temp_dipole.cpp @@ -15,6 +15,7 @@ #include "compute_temp_dipole.h" #include "atom.h" #include "force.h" +#include "domain.h" #include "group.h" #include "modify.h" #include "fix.h" @@ -76,7 +77,7 @@ void ComputeTempDipole::init() void ComputeTempDipole::recount() { double natoms = group->count(igroup); - dof = 2.0 * force->dimension * natoms; + dof = 2.0 * domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; diff --git a/src/GRANULAR/atom_vec_granular.cpp b/src/GRANULAR/atom_vec_granular.cpp index 744fc35968..29e63ee9ee 100644 --- a/src/GRANULAR/atom_vec_granular.cpp +++ b/src/GRANULAR/atom_vec_granular.cpp @@ -40,7 +40,7 @@ AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) : xcol_data = 5; atom->radius_flag = atom->density_flag = atom->rmass_flag = 1; - atom->xphi_flag = atom->omega_flag = atom->torque_flag = 1; + atom->xorient_flag = atom->omega_flag = atom->torque_flag = 1; PI = 4.0*atan(1.0); } @@ -76,8 +76,8 @@ void AtomVecGranular::grow(int n) rmass = atom->rmass = (double *) memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass"); - xphi = atom->xphi = - memory->grow_2d_double_array(atom->xphi,nmax,3,"atom:xphi"); + xorient = atom->xorient = + memory->grow_2d_double_array(atom->xorient,nmax,3,"atom:xorient"); omega = atom->omega = memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega"); torque = atom->torque = @@ -106,9 +106,9 @@ void AtomVecGranular::copy(int i, int j) radius[j] = radius[i]; density[j] = density[i]; rmass[j] = rmass[i]; - xphi[j][0] = xphi[i][0]; - xphi[j][1] = xphi[i][1]; - xphi[j][2] = xphi[i][2]; + xorient[j][0] = xorient[i][0]; + xorient[j][1] = xorient[i][1]; + xorient[j][2] = xorient[i][2]; omega[j][0] = omega[i][0]; omega[j][1] = omega[i][1]; omega[j][2] = omega[i][2]; @@ -407,9 +407,9 @@ int AtomVecGranular::pack_exchange(int i, double *buf) buf[m++] = radius[i]; buf[m++] = density[i]; buf[m++] = rmass[i]; - buf[m++] = xphi[i][0]; - buf[m++] = xphi[i][1]; - buf[m++] = xphi[i][2]; + buf[m++] = xorient[i][0]; + buf[m++] = xorient[i][1]; + buf[m++] = xorient[i][2]; buf[m++] = omega[i][0]; buf[m++] = omega[i][1]; buf[m++] = omega[i][2]; @@ -444,9 +444,9 @@ int AtomVecGranular::unpack_exchange(double *buf) radius[nlocal] = buf[m++]; density[nlocal] = buf[m++]; rmass[nlocal] = buf[m++]; - xphi[nlocal][0] = buf[m++]; - xphi[nlocal][1] = buf[m++]; - xphi[nlocal][2] = buf[m++]; + xorient[nlocal][0] = buf[m++]; + xorient[nlocal][1] = buf[m++]; + xorient[nlocal][2] = buf[m++]; omega[nlocal][0] = buf[m++]; omega[nlocal][1] = buf[m++]; omega[nlocal][2] = buf[m++]; @@ -502,9 +502,9 @@ int AtomVecGranular::pack_restart(int i, double *buf) buf[m++] = radius[i]; buf[m++] = density[i]; - buf[m++] = xphi[i][0]; - buf[m++] = xphi[i][1]; - buf[m++] = xphi[i][2]; + buf[m++] = xorient[i][0]; + buf[m++] = xorient[i][1]; + buf[m++] = xorient[i][2]; buf[m++] = omega[i][0]; buf[m++] = omega[i][1]; buf[m++] = omega[i][2]; @@ -546,15 +546,15 @@ int AtomVecGranular::unpack_restart(double *buf) radius[nlocal] = buf[m++]; density[nlocal] = buf[m++]; - if (force->dimension == 3) + if (domain->dimension == 3) rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; else rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; - xphi[nlocal][0] = buf[m++]; - xphi[nlocal][1] = buf[m++]; - xphi[nlocal][2] = buf[m++]; + xorient[nlocal][0] = buf[m++]; + xorient[nlocal][1] = buf[m++]; + xorient[nlocal][2] = buf[m++]; omega[nlocal][0] = buf[m++]; omega[nlocal][1] = buf[m++]; omega[nlocal][2] = buf[m++]; @@ -592,14 +592,14 @@ void AtomVecGranular::create_atom(int itype, double *coord) radius[nlocal] = 0.5; density[nlocal] = 1.0; - if (force->dimension == 3) + if (domain->dimension == 3) rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; else rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal]; - xphi[nlocal][0] = 0.0; - xphi[nlocal][1] = 0.0; - xphi[nlocal][2] = 0.0; + xorient[nlocal][0] = 0.0; + xorient[nlocal][1] = 0.0; + xorient[nlocal][2] = 0.0; omega[nlocal][0] = 0.0; omega[nlocal][1] = 0.0; omega[nlocal][2] = 0.0; @@ -627,7 +627,7 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) radius[nlocal] = 0.5 * atof(values[2]); density[nlocal] = atof(values[3]); - if (force->dimension == 3) + if (domain->dimension == 3) rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; else @@ -643,9 +643,9 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; - xphi[nlocal][0] = 0.0; - xphi[nlocal][1] = 0.0; - xphi[nlocal][2] = 0.0; + xorient[nlocal][0] = 1.0; + xorient[nlocal][1] = 0.0; + xorient[nlocal][2] = 0.0; omega[nlocal][0] = 0.0; omega[nlocal][1] = 0.0; omega[nlocal][2] = 0.0; @@ -662,7 +662,7 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) { radius[nlocal] = 0.5 * atof(values[0]); density[nlocal] = atof(values[1]); - if (force->dimension == 3) + if (domain->dimension == 3) rmass[nlocal] = 4.0*PI/3.0 * radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal]; else @@ -671,9 +671,9 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) v[nlocal][0] = 0.0; v[nlocal][1] = 0.0; v[nlocal][2] = 0.0; - xphi[nlocal][0] = 0.0; - xphi[nlocal][1] = 0.0; - xphi[nlocal][2] = 0.0; + xorient[nlocal][0] = 1.0; + xorient[nlocal][1] = 0.0; + xorient[nlocal][2] = 0.0; omega[nlocal][0] = 0.0; omega[nlocal][1] = 0.0; omega[nlocal][2] = 0.0; @@ -726,7 +726,7 @@ int AtomVecGranular::memory_usage() if (atom->memcheck("radius")) bytes += nmax * sizeof(double); if (atom->memcheck("density")) bytes += nmax * sizeof(double); if (atom->memcheck("rmass")) bytes += nmax * sizeof(double); - if (atom->memcheck("xphi")) bytes += nmax*3 * sizeof(double); + if (atom->memcheck("xorient")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double); if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double); diff --git a/src/GRANULAR/atom_vec_granular.h b/src/GRANULAR/atom_vec_granular.h index fc87ddf60b..b30a9a4567 100644 --- a/src/GRANULAR/atom_vec_granular.h +++ b/src/GRANULAR/atom_vec_granular.h @@ -53,7 +53,7 @@ class AtomVecGranular : public AtomVec { int *tag,*type,*mask,*image; double **x,**v,**f; double *radius,*density,*rmass; - double **xphi,**omega,**torque; + double **xorient,**omega,**torque; }; } diff --git a/src/GRANULAR/fix_nve_gran.cpp b/src/GRANULAR/fix_nve_gran.cpp index 3b026fb333..f169e40809 100644 --- a/src/GRANULAR/fix_nve_gran.cpp +++ b/src/GRANULAR/fix_nve_gran.cpp @@ -17,6 +17,7 @@ #include "atom.h" #include "update.h" #include "force.h" +#include "domain.h" #include "error.h" using namespace LAMMPS_NS; @@ -33,9 +34,9 @@ FixNVEGran::FixNVEGran(LAMMPS *lmp, int narg, char **arg) : { if (narg != 3) error->all("Illegal fix nve/gran command"); - if (!atom->xphi_flag || !atom->omega_flag || !atom->torque_flag) + if (!atom->xorient_flag || !atom->omega_flag || !atom->torque_flag) error->all("Fix nve/gran requires atom attributes " - "xphi, omega, torque"); + "xorient, omega, torque"); } /* ---------------------------------------------------------------------- */ @@ -54,8 +55,8 @@ void FixNVEGran::init() { dtv = update->dt; dtf = 0.5 * update->dt * force->ftm2v; - if (force->dimension == 3) dtfphi = dtf / INERTIA3D; - else dtfphi = dtf / INERTIA2D; + if (domain->dimension == 3) dtfrotate = dtf / INERTIA3D; + else dtfrotate = dtf / INERTIA2D; } /* ---------------------------------------------------------------------- */ @@ -67,7 +68,7 @@ void FixNVEGran::initial_integrate() double **x = atom->x; double **v = atom->v; double **f = atom->f; - double **xphi = atom->xphi; + double **xorient = atom->xorient; double **omega = atom->omega; double **torque = atom->torque; double *rmass = atom->rmass; @@ -85,13 +86,13 @@ void FixNVEGran::initial_integrate() x[i][1] += dtv * v[i][1]; x[i][2] += dtv * v[i][2]; - dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]); + dtfm = dtfrotate / (radius[i]*radius[i]*rmass[i]); omega[i][0] += dtfm * torque[i][0]; omega[i][1] += dtfm * torque[i][1]; omega[i][2] += dtfm * torque[i][2]; - xphi[i][0] += dtv * omega[i][0]; - xphi[i][1] += dtv * omega[i][1]; - xphi[i][2] += dtv * omega[i][2]; + xorient[i][0] += dtv * omega[i][0]; + xorient[i][1] += dtv * omega[i][1]; + xorient[i][2] += dtv * omega[i][2]; } } } @@ -118,7 +119,7 @@ void FixNVEGran::final_integrate() v[i][1] += dtfm * f[i][1]; v[i][2] += dtfm * f[i][2]; - dtfm = dtfphi / (radius[i]*radius[i]*rmass[i]); + dtfm = dtfrotate / (radius[i]*radius[i]*rmass[i]); omega[i][0] += dtfm * torque[i][0]; omega[i][1] += dtfm * torque[i][1]; omega[i][2] += dtfm * torque[i][2]; diff --git a/src/GRANULAR/fix_nve_gran.h b/src/GRANULAR/fix_nve_gran.h index 0dc4eaaf71..bb6644bacd 100644 --- a/src/GRANULAR/fix_nve_gran.h +++ b/src/GRANULAR/fix_nve_gran.h @@ -27,7 +27,7 @@ class FixNVEGran : public Fix { void final_integrate(); private: - double dtv,dtf,dtfphi; + double dtv,dtf,dtfrotate; }; } diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index f48da3e644..5d150e0447 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -94,7 +94,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : rate = atof(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"vel") == 0) { - if (force->dimension == 3) { + if (domain->dimension == 3) { if (iarg+6 > narg) error->all("Illegal fix pour command"); vxlo = atof(arg[iarg+1]); vxhi = atof(arg[iarg+2]); @@ -150,7 +150,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : error->all("Insertion region extends outside simulation box"); } else error->all("Must use a block or cylinder region with fix pour"); - if (region_style == 2 && force->dimension == 2) + if (region_style == 2 && domain->dimension == 2) error->all("Must use a block region with fix pour for 2d simulations"); // random number generator, same for all procs @@ -171,7 +171,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : double v_relative,delta; double g = 1.0; - if (force->dimension == 3) { + if (domain->dimension == 3) { v_relative = vz - rate; delta = v_relative + sqrt(v_relative*v_relative + 2.0*g*(zhi-zlo)) / g; } else { @@ -192,7 +192,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // in 3d, insure dy >= 1, for quasi-2d simulations double volume,volume_one; - if (force->dimension == 3) { + if (domain->dimension == 3) { if (region_style == 1) { double dy = yhi - ylo; if (dy < 1.0) dy = 1.0; @@ -261,7 +261,7 @@ void FixPour::init() double ygrav = sin(degree2rad * theta) * sin(degree2rad * phi); double zgrav = cos(degree2rad * theta); - if (force->dimension == 3) { + if (domain->dimension == 3) { if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON || fabs(zgrav+1.0) > EPSILON) error->all("Gravity must point in -z to use with fix pour in 3d"); @@ -299,7 +299,7 @@ void FixPour::pre_exchange() // lo/hi current = z (or y) bounds of insertion region this timestep - if (force->dimension == 3) { + if (domain->dimension == 3) { lo_current = zlo + (update->ntimestep - nfirst) * update->dt * rate; hi_current = zhi + (update->ntimestep - nfirst) * update->dt * rate; } else { @@ -430,7 +430,7 @@ void FixPour::pre_exchange() coord[2] = xnear[i][2]; radtmp = xnear[i][3]; denstmp = density_lo + random->uniform() * (density_hi-density_lo); - if (force->dimension == 3) { + if (domain->dimension == 3) { vxtmp = vxlo + random->uniform() * (vxhi-vxlo); vytmp = vylo + random->uniform() * (vyhi-vylo); vztmp = vz - sqrt(2.0*g*(hi_current-coord[2])); @@ -444,11 +444,11 @@ void FixPour::pre_exchange() if (coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1] && coord[2] >= sublo[2] && coord[2] < subhi[2]) flag = 1; - else if (force->dimension == 3 && coord[2] >= domain->boxhi[2] && + else if (domain->dimension == 3 && coord[2] >= domain->boxhi[2] && comm->myloc[2] == comm->procgrid[2]-1 && coord[0] >= sublo[0] && coord[0] < subhi[0] && coord[1] >= sublo[1] && coord[1] < subhi[1]) flag = 1; - else if (force->dimension == 2 && coord[1] >= domain->boxhi[1] && + else if (domain->dimension == 2 && coord[1] >= domain->boxhi[1] && comm->myloc[1] == comm->procgrid[1]-1 && coord[0] >= sublo[0] && coord[0] < subhi[0]) flag = 1; @@ -458,7 +458,7 @@ void FixPour::pre_exchange() atom->type[m] = ntype; atom->radius[m] = radtmp; atom->density[m] = denstmp; - if (force->dimension == 3) + if (domain->dimension == 3) atom->rmass[m] = 4.0*PI/3.0 * radtmp*radtmp*radtmp * denstmp; else atom->rmass[m] = PI * radtmp*radtmp * denstmp; @@ -505,7 +505,7 @@ int FixPour::overlap(int i) double delta = radius_hi + atom->radius[i]; double **x = atom->x; - if (force->dimension == 3) { + if (domain->dimension == 3) { if (region_style == 1) { if (x[i][0] < xlo-delta || x[i][0] > xhi+delta || x[i][1] < ylo-delta || x[i][1] > yhi+delta || @@ -530,7 +530,7 @@ int FixPour::overlap(int i) void FixPour::xyz_random(double h, double *coord) { - if (force->dimension == 3) { + if (domain->dimension == 3) { if (region_style == 1) { coord[0] = xlo + random->uniform() * (xhi-xlo); coord[1] = ylo + random->uniform() * (yhi-ylo); diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index c9dd68ed2a..c634460f9b 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -81,7 +81,7 @@ void Ewald::init() // error check if (domain->triclinic) error->all("Cannot use Ewald with triclinic box"); - if (force->dimension == 2) error->all("Cannot use Ewald with 2d simulation"); + if (domain->dimension == 2) error->all("Cannot use Ewald with 2d simulation"); if (!atom->q_flag) error->all("Kspace style requires atom attribute q"); diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index a687a17e79..e9cbe18f4c 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -107,7 +107,7 @@ void PPPM::init() if (domain->triclinic) error->all("Cannot (yet) use PPPM with triclinic box"); - if (force->dimension == 2) error->all("Cannot use PPPM with 2d simulation"); + if (domain->dimension == 2) error->all("Cannot use PPPM with 2d simulation"); if (!atom->q_flag) error->all("Kspace style requires atom attribute q"); diff --git a/src/atom.cpp b/src/atom.cpp index 2b02a844ee..97cf66558e 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -61,7 +61,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) molecule = NULL; q = NULL; mu = NULL; - xphi = quat = omega = angmom = torque = NULL; + xorient = quat = omega = angmom = torque = NULL; radius = density = rmass = vfrac = NULL; maxspecial = 1; @@ -87,7 +87,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) molecule_flag = 0; q_flag = mu_flag = 0; - xphi_flag = quat_flag = omega_flag = angmom_flag = torque_flag = 0; + xorient_flag = quat_flag = omega_flag = angmom_flag = torque_flag = 0; radius_flag = density_flag = rmass_flag = vfrac_flag = 0; // ntype-length arrays @@ -150,7 +150,7 @@ Atom::~Atom() memory->sfree(q); memory->destroy_2d_double_array(mu); - memory->destroy_2d_double_array(xphi); + memory->destroy_2d_double_array(xorient); memory->destroy_2d_double_array(quat); memory->destroy_2d_double_array(omega); memory->destroy_2d_double_array(angmom); diff --git a/src/atom.h b/src/atom.h index e46c024dcd..3f18c33c14 100644 --- a/src/atom.h +++ b/src/atom.h @@ -43,7 +43,7 @@ class Atom : protected Pointers { int *molecule; double *q,**mu; - double **xphi,**quat,**omega,**angmom,**torque; + double **xorient,**quat,**omega,**angmom,**torque; double *radius,*density,*rmass,*vfrac; int maxspecial; @@ -70,7 +70,7 @@ class Atom : protected Pointers { int molecule_flag; int q_flag,mu_flag; - int xphi_flag,quat_flag,omega_flag,angmom_flag,torque_flag; + int xorient_flag,quat_flag,omega_flag,angmom_flag,torque_flag; int radius_flag,density_flag,rmass_flag,vfrac_flag; // extra peratom info in restart file destined for fix & diag diff --git a/src/comm.cpp b/src/comm.cpp index c955e0ce14..efb58bb87d 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -166,7 +166,7 @@ void Comm::set_procs() if (procgrid[0]*procgrid[1]*procgrid[2] != nprocs) error->all("Bad grid of processors"); - if (force->dimension == 2 && procgrid[2] != 1) + if (domain->dimension == 2 && procgrid[2] != 1) error->all("Proc grid in z != 1 for 2d simulation"); if (grid2proc) memory->destroy_3d_int_array(grid2proc); @@ -272,7 +272,7 @@ void Comm::setup() need[0] = static_cast (cutghost[0] * procgrid[0] / prd[0]) + 1; need[1] = static_cast (cutghost[1] * procgrid[1] / prd[1]) + 1; need[2] = static_cast (cutghost[2] * procgrid[2] / prd[2]) + 1; - if (force->dimension == 2) need[2] = 0; + if (domain->dimension == 2) need[2] = 0; // if non-periodic, do not communicate further than procgrid-1 away // this enables very large cutoffs in non-periodic systems @@ -1300,7 +1300,7 @@ void Comm::procs2box() while (ipy <= nremain) { if (nremain % ipy == 0) { ipz = nremain/ipy; - if (force->dimension == 3 || ipz == 1) { + if (domain->dimension == 3 || ipz == 1) { surf = area[0]/ipx/ipy + area[1]/ipx/ipz + area[2]/ipy/ipz; if (surf < bestsurf) { bestsurf = surf; diff --git a/src/compute_rotate_dipole.cpp b/src/compute_rotate_dipole.cpp index f6bc4e4015..2830291aac 100644 --- a/src/compute_rotate_dipole.cpp +++ b/src/compute_rotate_dipole.cpp @@ -14,7 +14,7 @@ #include "mpi.h" #include "compute_rotate_dipole.h" #include "atom.h" -#include "force.h" +#include "domain.h" #include "group.h" #include "error.h" @@ -56,7 +56,7 @@ void ComputeRotateDipole::init() double *mass = atom->mass; double **shape = atom->shape; - if (force->dimension == 3) + if (domain->dimension == 3) for (int i = 1; i <= atom->ntypes; i++) inertia[i] = INERTIA3D * mass[i] * 0.25*shape[i][0]*shape[i][0]; else diff --git a/src/compute_rotate_gran.cpp b/src/compute_rotate_gran.cpp index 3f1c475770..a1cd49ef6f 100644 --- a/src/compute_rotate_gran.cpp +++ b/src/compute_rotate_gran.cpp @@ -15,6 +15,7 @@ #include "compute_rotate_gran.h" #include "atom.h" #include "force.h" +#include "domain.h" #include "group.h" #include "error.h" @@ -42,7 +43,7 @@ ComputeRotateGran::ComputeRotateGran(LAMMPS *lmp, int narg, char **arg) : void ComputeRotateGran::init() { - if (force->dimension == 3) pfactor = 0.5 * force->mvv2e * INERTIA3D; + if (domain->dimension == 3) pfactor = 0.5 * force->mvv2e * INERTIA3D; else pfactor = 0.5 * force->mvv2e * INERTIA2D; } diff --git a/src/compute_temp.cpp b/src/compute_temp.cpp index b4e2b2e6be..0806761bc2 100644 --- a/src/compute_temp.cpp +++ b/src/compute_temp.cpp @@ -15,6 +15,7 @@ #include "compute_temp.h" #include "atom.h" #include "force.h" +#include "domain.h" #include "modify.h" #include "fix.h" #include "group.h" @@ -59,7 +60,7 @@ void ComputeTemp::init() void ComputeTemp::recount() { double natoms = group->count(igroup); - dof = force->dimension * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; diff --git a/src/compute_temp_deform.cpp b/src/compute_temp_deform.cpp index 8c883d46db..bc48451155 100644 --- a/src/compute_temp_deform.cpp +++ b/src/compute_temp_deform.cpp @@ -83,7 +83,7 @@ void ComputeTempDeform::init() void ComputeTempDeform::recount() { double natoms = group->count(igroup); - dof = force->dimension * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; diff --git a/src/compute_temp_ramp.cpp b/src/compute_temp_ramp.cpp index d7025caa53..f0bb271ebb 100644 --- a/src/compute_temp_ramp.cpp +++ b/src/compute_temp_ramp.cpp @@ -129,7 +129,7 @@ void ComputeTempRamp::init() void ComputeTempRamp::recount() { double natoms = group->count(igroup); - dof = force->dimension * natoms; + dof = domain->dimension * natoms; dof -= extra_dof + fix_dof; if (dof > 0) tfactor = force->mvv2e / (dof * force->boltz); else tfactor = 0.0; diff --git a/src/compute_temp_region.cpp b/src/compute_temp_region.cpp index fa781859b0..0253187004 100644 --- a/src/compute_temp_region.cpp +++ b/src/compute_temp_region.cpp @@ -90,7 +90,7 @@ double ComputeTempRegion::compute_scalar() tarray[0] = count; tarray[1] = t; MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world); - dof = force->dimension * tarray_all[0] - extra_dof; + dof = domain->dimension * tarray_all[0] - extra_dof; if (dof > 0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); else scalar = 0.0; return scalar; diff --git a/src/create_box.cpp b/src/create_box.cpp index 398d28f968..d3778be6a5 100644 --- a/src/create_box.cpp +++ b/src/create_box.cpp @@ -37,7 +37,7 @@ void CreateBox::command(int narg, char **arg) if (domain->box_exist) error->all("Cannot create_box after simulation box is defined"); - if (force->dimension == 2 && domain->zperiodic == 0) + if (domain->dimension == 2 && domain->zperiodic == 0) error->all("Cannot run 2d simulation with nonperiodic Z dimension"); domain->box_exist = 1; diff --git a/src/domain.cpp b/src/domain.cpp index 2d7bd34fef..ce992bc16c 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -55,6 +55,7 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp) { box_exist = 0; + dimension = 3; nonperiodic = 0; xperiodic = yperiodic = zperiodic = 1; periodicity[0] = xperiodic; @@ -132,7 +133,7 @@ void Domain::set_initial_box() error->one("Box bounds are invalid"); if (triclinic) { - if (force->dimension == 2 && (xz != 0.0 || yz != 0.0)) + if (domain->dimension == 2 && (xz != 0.0 || yz != 0.0)) error->all("Cannot skew triclinic box in z for 2d simulation"); if (xy != 0.0 && (!xperiodic || !yperiodic)) error->all("Triclinic box must be periodic in skewed dimensions"); diff --git a/src/domain.h b/src/domain.h index 4f2a7c6005..4f31ef57db 100644 --- a/src/domain.h +++ b/src/domain.h @@ -21,7 +21,7 @@ namespace LAMMPS_NS { class Domain : protected Pointers { public: int box_exist; // 0 = not yet created, 1 = exists - + int dimension; // 2 = 2d, 3 = 3d int nonperiodic; // 0 = periodic in all 3 dims // 1 = periodic or fixed in all 6 // 2 = shrink-wrap in any of 6 diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 6cf71381d7..3456b3ab0d 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -570,7 +570,7 @@ void DumpCustom::parse_fields(int narg, char **arg) pack_choice[i] = &DumpCustom::pack_tag; vtype[i] = INT; } else if (strcmp(arg[iarg],"mol") == 0) { - if (atom->molecule == NULL) + if (!atom->molecule_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_molecule; vtype[i] = INT; @@ -641,58 +641,58 @@ void DumpCustom::parse_fields(int narg, char **arg) vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"q") == 0) { - if (atom->q == NULL) + if (!atom->q_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_q; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"mux") == 0) { - if (atom->mu == NULL) + if (!atom->mu_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_mux; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"muy") == 0) { - if (atom->mu == NULL) + if (!atom->mu_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"muz") == 0) { - if (atom->mu == NULL) + if (!atom->mu_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_muz; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quatw") == 0) { - if (atom->quat == NULL) + if (!atom->quat_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quatw; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quati") == 0) { - if (atom->quat == NULL) + if (!atom->quat_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quati; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quatj") == 0) { - if (atom->quat == NULL) + if (!atom->quat_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quatj; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"quatk") == 0) { - if (atom->quat == NULL) + if (!atom->quat_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_quatk; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqx") == 0) { - if (atom->torque == NULL) + if (!atom->torque_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqx; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqy") == 0) { - if (atom->torque == NULL) + if (!atom->torque_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqy; vtype[i] = DOUBLE; } else if (strcmp(arg[iarg],"tqz") == 0) { - if (atom->torque == NULL) + if (!atom->torque_flag) error->all("Dumping an atom quantity that isn't allocated"); pack_choice[i] = &DumpCustom::pack_tqz; vtype[i] = DOUBLE; diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 1e577a60d7..f78243fe05 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -162,20 +162,24 @@ FixDeform::FixDeform(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (set[5].style && (domain->xperiodic == 0 || domain->yperiodic == 0)) error->all("Cannot fix deform on a non-periodic boundary"); - // setup scaling + // apply scaling to FINAL,DELTA,VEL since they have distance/vel units - if (scaleflag && domain->lattice == NULL) + int flag = 0; + for (int i = 0; i < 6; i++) + if (set[i].style == FINAL || set[i].style == DELTA || + set[i].style == VEL) flag = 1; + + if (flag && scaleflag && domain->lattice == NULL) error->all("Use of fix deform with undefined lattice"); - if (scaleflag) { + if (flag && scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; zscale = domain->lattice->zlattice; } else xscale = yscale = zscale = 1.0; - // apply scaling to input parameters with dist/vel units - // for 3,4,5 scale is in 1st dimension, e.g. x for xz + // for 3,4,5 scaling is in 1st dimension, e.g. x for xz double map[6]; map[0] = xscale; map[1] = yscale; map[2] = zscale; diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp index 46a077e0a6..54b5393b3d 100644 --- a/src/fix_deposit.cpp +++ b/src/fix_deposit.cpp @@ -88,7 +88,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) : // apply scaling to all input parameters with dist/vel units - if (force->dimension == 2) { + if (domain->dimension == 2) { lo *= yscale; hi *= yscale; rate *= yscale; @@ -191,7 +191,7 @@ void FixDeposit::pre_exchange() // adjust vertical coord by offset - if (force->dimension == 2) coord[1] += offset; + if (domain->dimension == 2) coord[1] += offset; else coord[2] += offset; // if global, reset vertical coord to be lo-hi above highest atom @@ -202,7 +202,7 @@ void FixDeposit::pre_exchange() int dim; double max,maxall,delx,dely,delz,rsq; - if (force->dimension == 2) { + if (domain->dimension == 2) { dim = 1; max = domain->boxlo[1]; } else { @@ -218,7 +218,7 @@ void FixDeposit::pre_exchange() dely = coord[1] - x[i][1]; delz = 0.0; domain->minimum_image(delx,dely,delz); - if (force->dimension == 2) rsq = delx*delx; + if (domain->dimension == 2) rsq = delx*delx; else rsq = delx*delx + dely*dely; if (rsq > deltasq) continue; } @@ -226,7 +226,7 @@ void FixDeposit::pre_exchange() } MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world); - if (force->dimension == 2) + if (domain->dimension == 2) coord[1] = maxall + lo + random->uniform()*(hi-lo); else coord[2] = maxall + lo + random->uniform()*(hi-lo); @@ -271,11 +271,11 @@ void FixDeposit::pre_exchange() if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; - else if (force->dimension == 3 && newcoord[2] >= domain->boxhi[2] && + else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] && comm->myloc[2] == comm->procgrid[2]-1 && newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; - else if (force->dimension == 2 && newcoord[1] >= domain->boxhi[1] && + else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] && comm->myloc[1] == comm->procgrid[1]-1 && newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; diff --git a/src/fix_enforce2d.cpp b/src/fix_enforce2d.cpp index 90b4d91cff..2ca68ed17e 100644 --- a/src/fix_enforce2d.cpp +++ b/src/fix_enforce2d.cpp @@ -67,9 +67,7 @@ void FixEnforce2D::min_setup() void FixEnforce2D::post_force(int vflag) { double **v = atom->v; - double **omega = atom->omega; double **f = atom->f; - double **torque = atom->torque; int *mask = atom->mask; int nlocal = atom->nlocal; @@ -79,9 +77,10 @@ void FixEnforce2D::post_force(int vflag) f[i][2] = 0.0; } - // for omega/torque systems, zero xy rotation + // for systems with omega/angmom/torque, zero x and y components - if (omega) { + if (atom->omega_flag) { + double **omega = atom->omega; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { omega[i][0] = 0.0; @@ -89,7 +88,17 @@ void FixEnforce2D::post_force(int vflag) } } - if (torque) { + if (atom->angmom_flag) { + double **angmom = atom->angmom; + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + angmom[i][0] = 0.0; + angmom[i][1] = 0.0; + } + } + + if (atom->torque_flag) { + double **torque = atom->torque; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { torque[i][0] = 0.0; diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index 11aa45dbef..7768dbb137 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -71,7 +71,7 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) : newarg[1] = group->names[igroup]; if (strcmp(style,"nvt") == 0) newarg[2] = "temp"; else if (strcmp(style,"nvt/asphere") == 0) newarg[2] = "temp/asphere"; - else if (strcmp(style,"nvt/deform") == 0) newarg[2] = "temp/deform"; + else if (strcmp(style,"nvt/sllod") == 0) newarg[2] = "temp/deform"; modify->add_compute(3,newarg); delete [] newarg; tflag = 1; diff --git a/src/force.cpp b/src/force.cpp index 66b7a1ed14..802c1061d2 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -51,7 +51,6 @@ using namespace LAMMPS_NS; Force::Force(LAMMPS *lmp) : Pointers(lmp) { - dimension = 3; newton = newton_pair = newton_bond = 1; special_lj[1] = special_lj[2] = special_lj[3] = 0.0; diff --git a/src/force.h b/src/force.h index e3876d53eb..a1ae61ddd1 100644 --- a/src/force.h +++ b/src/force.h @@ -29,7 +29,6 @@ class Force : protected Pointers { double dielectric; // dielectric constant double qqrd2e; // q^2/r to energy w/ dielectric constant - int dimension; // 2 = 2d, 3 = 3d int newton,newton_pair,newton_bond; // Newton's 3rd law settings class Pair *pair; diff --git a/src/input.cpp b/src/input.cpp index 76f4387d7f..6a635fffb9 100644 --- a/src/input.cpp +++ b/src/input.cpp @@ -795,7 +795,7 @@ void Input::dimension() if (narg != 1) error->all("Illegal dimension command"); if (domain->box_exist) error->all("Dimension command after simulation box is defined"); - force->dimension = atoi(arg[0]); + domain->dimension = atoi(arg[0]); } /* ---------------------------------------------------------------------- */ diff --git a/src/lattice.cpp b/src/lattice.cpp index 9be04bad73..fe9c4f9b01 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -16,7 +16,7 @@ #include "stdlib.h" #include "lattice.h" #include "update.h" -#include "force.h" +#include "domain.h" #include "comm.h" #include "memory.h" #include "error.h" @@ -59,7 +59,7 @@ Lattice::Lattice(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) // check that lattice matches dimension // style CUSTOM can be either 2d or 3d - int dimension = force->dimension; + int dimension = domain->dimension; if (dimension == 2) { if (style == SC || style == BCC || style == FCC || style == DIAMOND) error->all("Lattice style incompatible with simulation dimension"); diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 21f93eb861..021574a6cd 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -520,19 +520,19 @@ void Neighbor::init() } else if (style == BIN) { if (force->newton_pair == 0) { half_build = &Neighbor::granular_bin_no_newton; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_no_newton; else half_stencil = &Neighbor::stencil_half_2d_no_newton; } else if (triclinic) { half_build = &Neighbor::granular_bin_newton_tri; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton_tri; else half_stencil = &Neighbor::stencil_half_2d_newton_tri; } else { half_build = &Neighbor::granular_bin_newton; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton; else half_stencil = &Neighbor::stencil_half_2d_newton; @@ -547,19 +547,19 @@ void Neighbor::init() } else if (style == BIN) { if (force->newton_pair == 0) { half_build = &Neighbor::respa_bin_no_newton; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_no_newton; else half_stencil = &Neighbor::stencil_half_2d_no_newton; } else if (triclinic) { half_build = &Neighbor::respa_bin_newton_tri; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton_tri; else half_stencil = &Neighbor::stencil_half_2d_newton_tri; } else { half_build = &Neighbor::respa_bin_newton; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton; else half_stencil = &Neighbor::stencil_half_2d_newton; @@ -582,7 +582,7 @@ void Neighbor::init() half_stencil = &Neighbor::stencil_none; } else { half_build = &Neighbor::half_bin_no_newton; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_no_newton; else half_stencil = &Neighbor::stencil_half_2d_no_newton; @@ -593,13 +593,13 @@ void Neighbor::init() half_stencil = &Neighbor::stencil_none; } else if (triclinic) { half_build = &Neighbor::half_bin_newton_tri; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton_tri; else half_stencil = &Neighbor::stencil_half_2d_newton_tri; } else { half_build = &Neighbor::half_bin_newton; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton; else half_stencil = &Neighbor::stencil_half_2d_newton; @@ -612,7 +612,7 @@ void Neighbor::init() half_stencil = &Neighbor::stencil_none; } else { half_build = &Neighbor::half_bin_no_newton_multi; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_no_newton_multi; else half_stencil = &Neighbor::stencil_half_2d_no_newton_multi; @@ -623,13 +623,13 @@ void Neighbor::init() half_stencil = &Neighbor::stencil_none; } else if (triclinic) { half_build = &Neighbor::half_bin_newton_multi_tri; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton_multi_tri; else half_stencil = &Neighbor::stencil_half_2d_newton_multi_tri; } else { half_build = &Neighbor::half_bin_newton_multi; - if (force->dimension == 3) + if (domain->dimension == 3) half_stencil = &Neighbor::stencil_half_3d_newton_multi; else half_stencil = &Neighbor::stencil_half_2d_newton_multi; @@ -644,13 +644,13 @@ void Neighbor::init() if (style == NSQ) full_build = &Neighbor::full_nsq; else if (style == BIN) { full_build = &Neighbor::full_bin; - if (force->dimension == 3) + if (domain->dimension == 3) full_stencil = &Neighbor::stencil_full_3d; else full_stencil = &Neighbor::stencil_full_2d; } else { full_build = &Neighbor::full_bin_multi; - if (force->dimension == 3) + if (domain->dimension == 3) full_stencil = &Neighbor::stencil_full_3d_multi; else full_stencil = &Neighbor::stencil_full_2d_multi; @@ -1001,7 +1001,7 @@ void Neighbor::setup_bins() nbinx = static_cast (bbox[0]*binsizeinv); nbiny = static_cast (bbox[1]*binsizeinv); - if (force->dimension == 3) + if (domain->dimension == 3) nbinz = static_cast (bbox[2]*binsizeinv); else nbinz = 1; diff --git a/src/pair.cpp b/src/pair.cpp index 1bbadd920c..98ca74e888 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -117,7 +117,7 @@ void Pair::init() if (offset_flag && tail_flag) error->all("Cannot have both pair_modify shift and tail set to yes"); - if (tail_flag && force->dimension == 2) + if (tail_flag && domain->dimension == 2) error->all("Cannot use pair tail corrections with 2d simulations"); if (tail_flag && domain->nonperiodic && comm->me == 0) error->warning("Using pair tail corrections with nonperiodic system"); diff --git a/src/read_data.cpp b/src/read_data.cpp index 40e4dfdc18..b812253c7c 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -71,7 +71,7 @@ void ReadData::command(int narg, char **arg) if (domain->box_exist) error->all("Cannot read_data after simulation box is defined"); - if (force->dimension == 2 && domain->zperiodic == 0) + if (domain->dimension == 2 && domain->zperiodic == 0) error->all("Cannot run 2d simulation with nonperiodic Z dimension"); // scan data file to determine max topology needed per atom diff --git a/src/read_restart.cpp b/src/read_restart.cpp index b64ba974bb..2182f96937 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -378,10 +378,10 @@ void ReadRestart::header() } else if (flag == DIMENSION) { int dimension = read_int(); - if (dimension != force->dimension && me == 0) + if (dimension != domain->dimension && me == 0) error->warning("Resetting dimension to restart file value"); - force->dimension = dimension; - if (force->dimension == 2 && domain->zperiodic == 0) + domain->dimension = dimension; + if (domain->dimension == 2 && domain->zperiodic == 0) error->all("Cannot run 2d simulation with nonperiodic Z dimension"); // read nprocs_file from restart file, warn if different diff --git a/src/replicate.cpp b/src/replicate.cpp index 9b0fe26b94..401d994f9e 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -62,7 +62,7 @@ void Replicate::command(int narg, char **arg) // error and warning checks if (nx <= 0 || ny <= 0 || nz <= 0) error->all("Illegal replicate command"); - if (force->dimension == 2 && nz != 1) + if (domain->dimension == 2 && nz != 1) error->all("Cannot replicate 2d simulation in z dimension"); if ((nx > 1 && domain->xperiodic == 0) || (ny > 1 && domain->yperiodic == 0) || diff --git a/src/set.cpp b/src/set.cpp index 617d9c1510..3cea4cee24 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -26,6 +26,7 @@ #include "force.h" #include "pair.h" #include "random_park.h" +#include "math_extra.h" #include "error.h" using namespace LAMMPS_NS; @@ -94,7 +95,7 @@ void Set::command(int narg, char **arg) } else if (strcmp(arg[iarg],"mol") == 0) { if (iarg+2 > narg) error->all("Illegal set command"); ivalue = atoi(arg[iarg+1]); - if (atom->molecule == NULL) + if (!atom->molecule_flag) error->all("Cannot set this attribute for this atom style"); set(MOLECULE); iarg += 2; @@ -131,7 +132,7 @@ void Set::command(int narg, char **arg) } else if (strcmp(arg[iarg],"charge") == 0) { if (iarg+2 > narg) error->all("Illegal set command"); dvalue = atof(arg[iarg+1]); - if (atom->q == NULL) + if (!atom->q_flag) error->all("Cannot set this attribute for this atom style"); set(CHARGE); iarg += 2; @@ -140,31 +141,32 @@ void Set::command(int narg, char **arg) xvalue = atof(arg[iarg+1]); yvalue = atof(arg[iarg+2]); zvalue = atof(arg[iarg+3]); - if (atom->mu == NULL || atom->dipole == NULL) + if (!atom->mu_flag || atom->dipole == NULL) error->all("Cannot set this attribute for this atom style"); set(DIPOLE); iarg += 4; } else if (strcmp(arg[iarg],"dipole/random") == 0) { if (iarg+2 > narg) error->all("Illegal set command"); ivalue = atoi(arg[iarg+1]); - if (atom->mu == NULL || atom->shape == NULL) + if (!atom->mu_flag || atom->shape == NULL) error->all("Cannot set this attribute for this atom style"); if (ivalue <= 0) error->all("Invalid random number seed in set command"); setrandom(DIPOLE_RANDOM); iarg += 2; } else if (strcmp(arg[iarg],"quat") == 0) { - if (iarg+4 > narg) error->all("Illegal set command"); + if (iarg+5 > narg) error->all("Illegal set command"); xvalue = atof(arg[iarg+1]); yvalue = atof(arg[iarg+2]); zvalue = atof(arg[iarg+3]); - if (atom->quat == NULL) + wvalue = atof(arg[iarg+4]); + if (!atom->quat_flag) error->all("Cannot set this attribute for this atom style"); set(QUAT); - iarg += 4; + iarg += 5; } else if (strcmp(arg[iarg],"quat/random") == 0) { if (iarg+2 > narg) error->all("Illegal set command"); ivalue = atoi(arg[iarg+1]); - if (atom->quat == NULL) + if (!atom->quat_flag) error->all("Cannot set this attribute for this atom style"); if (ivalue <= 0) error->all("Invalid random number seed in set command"); setrandom(QUAT_RANDOM); @@ -305,11 +307,15 @@ void Set::set(int keyword) } } else if (keyword == QUAT) { - // need to set orientation correctly and normalize quat + double PI = 4.0*atan(1.0); + double theta2 = 0.5 * PI * wvalue/180.0; + double sintheta2 = sin(theta2); double **quat = atom->quat; - quat[i][0] = xvalue; - quat[i][1] = yvalue; - quat[i][2] = zvalue; + quat[i][0] = cos(theta2); + quat[i][1] = xvalue * sintheta2; + quat[i][2] = yvalue * sintheta2; + quat[i][3] = zvalue * sintheta2; + MathExtra::normalize4(quat[i]); } count++; } @@ -345,55 +351,90 @@ void Set::setrandom(int keyword) count++; } - // set dipole moments to random orientations + // set dipole moments to random orientations in 3d or 2d // dipole length is determined by dipole type array - } else if (keyword == DIPOLE) { + } else if (keyword == DIPOLE_RANDOM) { int *type = atom->type; double *dipole = atom->dipole; double **mu = atom->mu; int nlocal = atom->nlocal; double msq,scale; - for (i = 0; i < nlocal; i++) - if (select[i]) { - if (dipole[type[i]] > 0.0) { - random->reset(seed,x[i]); - mu[i][0] = random->uniform() - 0.5; - mu[i][1] = random->uniform() - 0.5; - mu[i][2] = random->uniform() - 0.5; - msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; - scale = dipole[type[i]]/sqrt(msq); - mu[i][0] *= scale; - mu[i][1] *= scale; - mu[i][2] *= scale; - count++; + + if (domain->dimension == 3) { + for (i = 0; i < nlocal; i++) + if (select[i]) { + if (dipole[type[i]] > 0.0) { + random->reset(seed,x[i]); + mu[i][0] = random->uniform() - 0.5; + mu[i][1] = random->uniform() - 0.5; + mu[i][2] = random->uniform() - 0.5; + msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; + scale = dipole[type[i]]/sqrt(msq); + mu[i][0] *= scale; + mu[i][1] *= scale; + mu[i][2] *= scale; + count++; + } } - } - // set quaternions to random orientations + } else { + for (i = 0; i < nlocal; i++) + if (select[i]) { + if (dipole[type[i]] > 0.0) { + random->reset(seed,x[i]); + mu[i][0] = random->uniform() - 0.5; + mu[i][1] = random->uniform() - 0.5; + mu[i][2] = 0.0; + msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1]; + scale = dipole[type[i]]/sqrt(msq); + mu[i][0] *= scale; + mu[i][1] *= scale; + count++; + } + } + } - } else if (keyword == QUAT) { + // set quaternions to random orientations in 3d or 2d + // no need to normalize quats since creations algorithms already do + + } else if (keyword == QUAT_RANDOM) { double **quat = atom->quat; int nlocal = atom->nlocal; - double s,t1,t2,theta1,theta2; - double PI = 4.0*atan(1.0); + if (domain->dimension == 3) { + double s,t1,t2,theta1,theta2; + double PI = 4.0*atan(1.0); + for (i = 0; i < nlocal; i++) + if (select[i]) { + random->reset(seed,x[i]); + s = random->uniform(); + t1 = sqrt(1.0-s); + t2 = sqrt(s); + theta1 = 2.0*PI*random->uniform(); + theta2 = 2.0*PI*random->uniform(); + quat[i][0] = cos(theta2)*t2; + quat[i][1] = sin(theta1)*t1; + quat[i][2] = cos(theta1)*t1; + quat[i][3] = sin(theta2)*t2; + count++; + } - for (i = 0; i < nlocal; i++) - if (select[i]) { - random->reset(seed,x[i]); - s = random->uniform(); - t1 = sqrt(1.0-s); - t2 = sqrt(s); - theta1 = 2.0*PI*random->uniform(); - theta2 = 2.0*PI*random->uniform(); - quat[i][0] = cos(theta2)*t2; - quat[i][1] = sin(theta1)*t1; - quat[i][2] = cos(theta1)*t1; - quat[i][3] = sin(theta2)*t2; - count++; - } + } else { + double theta2; + double PI = 4.0*atan(1.0); + for (i = 0; i < nlocal; i++) + if (select[i]) { + random->reset(seed,x[i]); + theta2 = PI*random->uniform(); + quat[i][0] = cos(theta2); + quat[i][1] = 0.0; + quat[i][2] = 0.0; + quat[i][3] = sin(theta2); + count++; + } + } } delete random; diff --git a/src/set.h b/src/set.h index f5327b117f..1136b8b0b7 100644 --- a/src/set.h +++ b/src/set.h @@ -27,7 +27,7 @@ class Set : protected Pointers { char *id; int *select; int style,ivalue,newtype,count; - double dvalue,xvalue,yvalue,zvalue,fraction; + double dvalue,xvalue,yvalue,zvalue,wvalue,fraction; void selection(int); void set(int); diff --git a/src/style.h b/src/style.h index da988dd2bb..8508a60a17 100644 --- a/src/style.h +++ b/src/style.h @@ -154,6 +154,7 @@ DumpStyle(xyz,DumpXYZ) #include "fix_nve.h" #include "fix_nve_noforce.h" #include "fix_nvt.h" +#include "fix_nvt_sllod.h" #include "fix_plane_force.h" #include "fix_print.h" #include "fix_orient_fcc.h" @@ -201,6 +202,7 @@ FixStyle(npt,FixNPT) FixStyle(nve,FixNVE) FixStyle(nve,FixNVENoforce) FixStyle(nvt,FixNVT) +FixStyle(nvt/sllod,FixNVTSlodd) FixStyle(orient/fcc,FixOrientFCC) FixStyle(print,FixPrint) FixStyle(planeforce,FixPlaneForce) diff --git a/src/velocity.cpp b/src/velocity.cpp index ab400ad78d..0ccf28c4b0 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -176,7 +176,7 @@ void Velocity::create(int narg, char **arg) int *mask = atom->mask; double *mass = atom->mass; double *rmass = atom->rmass; - int dimension = force->dimension; + int dimension = domain->dimension; int m; double vx,vy,vz,factor; @@ -360,7 +360,7 @@ void Velocity::set(int narg, char **arg) double **v = atom->v; int *mask = atom->mask; int nlocal = atom->nlocal; - int dimension = force->dimension; + int dimension = domain->dimension; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -427,7 +427,7 @@ void Velocity::ramp(int narg, char **arg) else if (strcmp(arg[0],"vz") == 0) v_dim = 2; else error->all("Illegal velocity command"); - if (v_dim == 2 && force->dimension == 2) + if (v_dim == 2 && domain->dimension == 2) error->all("Velocity ramp in z for a 2d problem"); double v_lo,v_hi; diff --git a/src/write_restart.cpp b/src/write_restart.cpp index 9c67385369..05fa1ac7ef 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -242,7 +242,7 @@ void WriteRestart::header() write_char(VERSION,universe->version); write_char(UNITS,update->unit_style); write_int(NTIMESTEP,update->ntimestep); - write_int(DIMENSION,force->dimension); + write_int(DIMENSION,domain->dimension); write_int(NPROCS,nprocs); write_int(PROCGRID_0,comm->procgrid[0]); write_int(PROCGRID_1,comm->procgrid[1]);