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@642 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2007-06-22 16:59:17 +00:00
parent 988a4479e8
commit cea0208207
42 changed files with 245 additions and 194 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/ASPHERE/atom_vec_ellipsoid.cpp
View File

@ -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

1
src/ASPHERE/atom_vec_ellipsoid.h
View File

@ -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 *);

3
src/ASPHERE/compute_temp_asphere.cpp
View File

@ -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

3
src/DIPOLE/compute_temp_dipole.cpp
View File

@ -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;

64
src/GRANULAR/atom_vec_granular.cpp
View File

@ -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);

2
src/GRANULAR/atom_vec_granular.h
View File

@ -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;
};
}

21
src/GRANULAR/fix_nve_gran.cpp
View File

@ -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];

2
src/GRANULAR/fix_nve_gran.h
View File

@ -27,7 +27,7 @@ class FixNVEGran : public Fix {
void final_integrate();
private:
double dtv,dtf,dtfphi;
double dtv,dtf,dtfrotate;
};
}

24
src/GRANULAR/fix_pour.cpp
View File

@ -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);

2
src/KSPACE/ewald.cpp
View File

@ -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");

2
src/KSPACE/pppm.cpp
View File

@ -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");

6
src/atom.cpp
View File

@ -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);

4
src/atom.h
View File

@ -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

6
src/comm.cpp
View File

@ -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<int> (cutghost[0] * procgrid[0] / prd[0]) + 1;
need[1] = static_cast<int> (cutghost[1] * procgrid[1] / prd[1]) + 1;
need[2] = static_cast<int> (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;

4
src/compute_rotate_dipole.cpp
View File

@ -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

3
src/compute_rotate_gran.cpp
View File

@ -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;
}

3
src/compute_temp.cpp
View File

@ -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;

2
src/compute_temp_deform.cpp
View File

@ -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;

2
src/compute_temp_ramp.cpp
View File

@ -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;

2
src/compute_temp_region.cpp
View File

@ -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;

2
src/create_box.cpp
View File

@ -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;

3
src/domain.cpp
View File

@ -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");

2
src/domain.h
View File

@ -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

24
src/dump_custom.cpp
View File

@ -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;

14
src/fix_deform.cpp
View File

@ -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;

14
src/fix_deposit.cpp
View File

@ -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;

19
src/fix_enforce2d.cpp
View File

@ -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;

2
src/fix_nvt.cpp
View File

@ -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;

1
src/force.cpp
View File

@ -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;

1
src/force.h
View File

@ -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;

2
src/input.cpp
View File

@ -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]);
}
/* ---------------------------------------------------------------------- */

4
src/lattice.cpp
View File

@ -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");

30
src/neighbor.cpp
View File

@ -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<int> (bbox[0]*binsizeinv);
nbiny = static_cast<int> (bbox[1]*binsizeinv);
if (force->dimension == 3)
if (domain->dimension == 3)
nbinz = static_cast<int> (bbox[2]*binsizeinv);
else nbinz = 1;

2
src/pair.cpp
View File

@ -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");

2
src/read_data.cpp
View File

@ -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

6
src/read_restart.cpp
View File

@ -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

2
src/replicate.cpp
View File

@ -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) ||

75
src/set.cpp
View File

@ -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,16 +351,18 @@ 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;
if (domain->dimension == 3) {
for (i = 0; i < nlocal; i++)
if (select[i]) {
if (dipole[type[i]] > 0.0) {
@ -371,15 +379,33 @@ void Set::setrandom(int keyword)
}
}
// 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;
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]);
@ -394,6 +420,21 @@ void Set::setrandom(int keyword)
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;

2
src/set.h
View File

@ -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);

2
src/style.h
View File

@ -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)

6
src/velocity.cpp
View File

@ -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;

2
src/write_restart.cpp
View File

@ -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]);