From c503b8736d386bbcff1133cfb232480a35f36107 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 3 Sep 2015 06:51:55 -0400 Subject: [PATCH] undo cosmetic differences between LAMMPS-ICMS and LAMMPS --- src/MPIIO/restart_mpiio.cpp | 1 - src/REPLICA/verlet_split.cpp | 2 +- src/RIGID/fix_rigid_nh.h | 2 +- src/USER-MISC/compute_basal_atom.cpp | 5 +- src/USER-MISC/fix_ipi.cpp | 4 +- src/USER-MISC/fix_srp.cpp | 70 +++++----- src/USER-MISC/fix_srp.h | 2 +- src/USER-MISC/pair_srp.cpp | 184 +++++++++++++-------------- src/USER-MISC/pair_srp.h | 4 +- src/USER-OMP/fix_nh_asphere_omp.cpp | 4 +- src/USER-OMP/fix_nh_sphere_omp.cpp | 2 +- src/USER-OMP/fix_rigid_nh_omp.cpp | 100 +++++++-------- src/USER-OMP/fix_rigid_nh_omp.h | 4 +- src/USER-OMP/fix_rigid_nph_omp.cpp | 16 +-- src/USER-OMP/fix_rigid_nph_omp.h | 2 +- src/USER-OMP/fix_rigid_npt_omp.cpp | 16 +-- src/USER-OMP/fix_rigid_npt_omp.h | 2 +- src/USER-OMP/fix_rigid_nve_omp.cpp | 2 +- src/USER-OMP/fix_rigid_nve_omp.h | 2 +- src/USER-OMP/fix_rigid_nvt_omp.cpp | 10 +- src/USER-OMP/fix_rigid_nvt_omp.h | 2 +- src/USER-OMP/fix_rigid_omp.cpp | 8 +- src/USER-OMP/fix_rigid_small_omp.cpp | 10 +- src/USER-OMP/thr_omp.cpp | 1 - src/USER-PHONON/fix_phonon.cpp | 8 +- src/atom.cpp | 16 +-- src/atom_vec_body.cpp | 36 +++--- src/delete_bonds.cpp | 2 - src/error.h | 1 - src/fix_ave_spatial.cpp | 2 +- src/fix_heat.cpp | 2 +- src/fix_nh.h | 3 +- src/force.h | 2 +- src/image.cpp | 4 +- src/lattice.cpp | 2 +- src/memory.h | 1 - src/pair.cpp | 4 +- src/pair.h | 6 +- src/random_mars.cpp | 4 +- src/random_park.cpp | 4 +- src/universe.h | 2 +- src/update.h | 2 +- 42 files changed, 268 insertions(+), 288 deletions(-) diff --git a/src/MPIIO/restart_mpiio.cpp b/src/MPIIO/restart_mpiio.cpp index 7846f0909a..b769a58a8a 100644 --- a/src/MPIIO/restart_mpiio.cpp +++ b/src/MPIIO/restart_mpiio.cpp @@ -16,7 +16,6 @@ ------------------------------------------------------------------------- */ #include "mpi.h" -#include "mpiio.h" #include "restart_mpiio.h" #include "error.h" #include "limits.h" diff --git a/src/REPLICA/verlet_split.cpp b/src/REPLICA/verlet_split.cpp index 2c79f1d421..408821fe22 100644 --- a/src/REPLICA/verlet_split.cpp +++ b/src/REPLICA/verlet_split.cpp @@ -194,7 +194,7 @@ VerletSplit::VerletSplit(LAMMPS *lmp, int narg, char **arg) : // f_kspace = Rspace copy of Kspace forces // allocate dummy version for Kspace partition - maxatom = -1; + maxatom = 0; f_kspace = NULL; if (!master) memory->create(f_kspace,1,1,"verlet/split:f_kspace"); } diff --git a/src/RIGID/fix_rigid_nh.h b/src/RIGID/fix_rigid_nh.h index f562475d50..64bfa37f90 100644 --- a/src/RIGID/fix_rigid_nh.h +++ b/src/RIGID/fix_rigid_nh.h @@ -72,7 +72,7 @@ class FixRigidNH : public FixRigid { int tcomputeflag,pcomputeflag; void couple(); - virtual void remap(); + void remap(); void nhc_temp_integrate(); void nhc_press_integrate(); diff --git a/src/USER-MISC/compute_basal_atom.cpp b/src/USER-MISC/compute_basal_atom.cpp index 891541263f..85a08d1838 100644 --- a/src/USER-MISC/compute_basal_atom.cpp +++ b/src/USER-MISC/compute_basal_atom.cpp @@ -204,10 +204,8 @@ void ComputeBasalAtom::compute_peratom() double bond_angle; double norm_j, norm_k; chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0; - double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik, + double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik,x3[n0],y3[n0],z3[n0], xmean5, ymean5, zmean5, xmean6, ymean6, zmean6, xmean7, ymean7, zmean7; - double *x3,*y3,*z3; - x3 = new double[n0]; y3 = new double[n0]; z3 = new double[n0]; for (j = 0; j < n0; j++) { x_ij = x[i][0]-x[nearest_n0[j]][0]; y_ij = x[i][1]-x[nearest_n0[j]][1]; @@ -415,7 +413,6 @@ void ComputeBasalAtom::compute_peratom() //if there are less than two ~180 degree bond angles, the algorithm returns null else BPV[i][0] = BPV[i][1] = BPV[i][2] = 0.0; - delete[] x3; delete[] y3; delete[] z3; //normalize BPV: double Mag = sqrt(BPV[i][0]*BPV[i][0] + BPV[i][1]*BPV[i][1] + BPV[i][2]*BPV[i][2]); diff --git a/src/USER-MISC/fix_ipi.cpp b/src/USER-MISC/fix_ipi.cpp index 77387fa6ed..58251a814e 100644 --- a/src/USER-MISC/fix_ipi.cpp +++ b/src/USER-MISC/fix_ipi.cpp @@ -380,8 +380,7 @@ void FixIPI::final_integrate() int nat=bsize/3; double **f= atom->f; - double *lbuf; - lbuf = new double[bsize]; + double lbuf[bsize]; // reassembles the force vector from the local arrays int nlocal = atom->nlocal; @@ -393,7 +392,6 @@ void FixIPI::final_integrate() lbuf[3*(atom->tag[i]-1)+2]=f[i][2]*forceconv; } MPI_Allreduce(lbuf,buffer,bsize,MPI_DOUBLE,MPI_SUM,world); - delete[] lbuf; for (int i = 0; i < 9; ++i) vir[i]=0.0; diff --git a/src/USER-MISC/fix_srp.cpp b/src/USER-MISC/fix_srp.cpp index ced1e6ee5e..08224d2950 100644 --- a/src/USER-MISC/fix_srp.cpp +++ b/src/USER-MISC/fix_srp.cpp @@ -49,22 +49,22 @@ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // per-atom array width 2 peratom_flag = 1; - size_peratom_cols = 2; + size_peratom_cols = 2; - // initial allocation of atom-based array + // initial allocation of atom-based array // register with Atom class array = NULL; - grow_arrays(atom->nmax); + grow_arrays(atom->nmax); - // extends pack_exchange() + // extends pack_exchange() atom->add_callback(0); atom->add_callback(1); // restart - atom->add_callback(2); + atom->add_callback(2); - // zero - for (int i = 0; i < atom->nmax; i++) - for (int m = 0; m < 3; m++) - array[i][m] = 0.0; + // zero + for (int i = 0; i < atom->nmax; i++) + for (int m = 0; m < 3; m++) + array[i][m] = 0.0; // assume setup of fix is needed to insert particles // is true if reading from datafile @@ -72,7 +72,7 @@ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) // might be true if reading from restart // a restart file written during the run has bond particles as per atom data // a restart file written after the run does not have bond particles - setup = 1; + setup = 1; } /* ---------------------------------------------------------------------- */ @@ -106,7 +106,7 @@ void FixSRP::init() error->all(FLERR,"Cannot use pair srp without pair_style hybrid"); // the highest numbered atom type should be reserved for bond particles (bptype) - // set bptype, unless it will be read from restart + // set bptype, unless it will be read from restart if(!restart_reset) bptype = atom->ntypes; // check if bptype is already in use @@ -118,24 +118,24 @@ void FixSRP::init() if(!restart_reset) error->all(FLERR,"Fix SRP requires an extra atom type"); else - setup = 0; + setup = 0; } - + // setup neigh exclusions for diff atom types // bond particles do not interact with other types // type bptype only interacts with itself char* arg1[4]; - arg1[0] = (char *) "exclude"; - arg1[1] = (char *) "type"; + arg1[0] = (char *) "exclude"; + arg1[1] = (char *) "type"; char c0[20]; char c1[20]; for(int z = 1; z < bptype; z++) { sprintf(c0, "%d", z); - arg1[2] = c0; + arg1[2] = c0; sprintf(c1, "%d", bptype); - arg1[3] = c1; + arg1[3] = c1; neighbor->modify_params(4, arg1); } } @@ -192,19 +192,19 @@ void FixSRP::setup_pre_force(int zz) xone[1] = (xold[i][1] + xold[j][1])*0.5; xone[2] = (xold[i][2] + xold[j][2])*0.5; - // record longest bond + // record longest bond // this used to set ghost cutoff delx = xold[j][0] - xold[i][0]; dely = xold[j][1] - xold[i][1]; delz = xold[j][2] - xold[i][2]; rsq = delx*delx + dely*dely + delz*delz; - if(rsq > rsqold) rsqold = rsq; + if(rsq > rsqold) rsqold = rsq; // make one particle for each bond // i is local // if newton bond, always make particle // if j is local, always make particle - // if j is ghost, decide from tag + // if j is ghost, decide from tag if( force->newton_bond || j < nlocal || tagold[i] > tagold[j] ){ atom->natoms += 1; @@ -246,7 +246,7 @@ void FixSRP::setup_pre_force(int zz) // ghost atoms must be present for bonds on edge of neighbor cutoff // extend cutghost slightly more than half of the longest bond MPI_Allreduce(&rsqold,&rsqmax,1,MPI_DOUBLE,MPI_MAX,world); - rmax = sqrt(rsqmax); + rmax = sqrt(rsqmax); double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax; // find smallest cutghost @@ -262,7 +262,7 @@ void FixSRP::setup_pre_force(int zz) sprintf(str, "Extending ghost comm cutoff. New %f, old %f.", cutneighmax_srp, cutghostmin); error->message(FLERR,str); } - // cutghost updated by comm->setup + // cutghost updated by comm->setup comm->cutghostuser = cutneighmax_srp; } @@ -270,7 +270,7 @@ void FixSRP::setup_pre_force(int zz) // move owned to new procs // get ghosts // build neigh lists again - + // if triclinic, lambda coords needed for pbc, exchange, borders if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); @@ -290,7 +290,7 @@ void FixSRP::setup_pre_force(int zz) // new atom counts nlocal = atom->nlocal; nall = atom->nlocal + atom->nghost; - // zero all forces + // zero all forces for(int i = 0; i < nall; i++) for(int n = 0; n < 3; n++) atom->f[i][n] = 0.0; @@ -308,7 +308,7 @@ void FixSRP::setup_pre_force(int zz) void FixSRP::pre_exchange() { - // update ghosts + // update ghosts comm->forward_comm(); // reassign bond particle coordinates to midpoint of bonds @@ -328,7 +328,7 @@ void FixSRP::pre_exchange() if(j < 0) error->all(FLERR,"Fix SRP failed to map atom"); j = domain->closest_image(ii,j); - // position of bond particle ii + // position of bond particle ii atom->x[ii][0] = (x[i][0] + x[j][0])*0.5; atom->x[ii][1] = (x[i][1] + x[j][1])*0.5; atom->x[ii][2] = (x[i][2] + x[j][2])*0.5; @@ -374,7 +374,7 @@ void FixSRP::copy_arrays(int i, int j, int delflag) void FixSRP::set_arrays(int i) { array[i][0] = -1; - array[i][1] = -1; + array[i][1] = -1; } /* ---------------------------------------------------------------------- @@ -438,9 +438,9 @@ int FixSRP::unpack_border(int n, int first, double *buf) void FixSRP::post_run() { // all bond particles are removed after each run - // useful for write_data and write_restart commands + // useful for write_data and write_restart commands // since those commands occur between runs - + bigint natoms_previous = atom->natoms; int nlocal = atom->nlocal; int* dlist; @@ -508,7 +508,7 @@ void FixSRP::post_run() comm->borders(); // change back to box coordinates if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); -} +} /* ---------------------------------------------------------------------- pack values in local atom-based arrays for restart file @@ -518,7 +518,7 @@ int FixSRP::pack_restart(int i, double *buf) { int m = 0; buf[m++] = 3; - buf[m++] = array[i][0]; + buf[m++] = array[i][0]; buf[m++] = array[i][1]; return m; } @@ -527,17 +527,17 @@ int FixSRP::pack_restart(int i, double *buf) unpack values from atom->extra array to restart the fix ------------------------------------------------------------------------- */ -void FixSRP::unpack_restart(int nlocal, int nth) +void FixSRP::unpack_restart(int nlocal, int nth) { double **extra = atom->extra; // skip to Nth set of extra values - int m = 0; + int m = 0; for (int i = 0; i < nth; i++){ m += static_cast (extra[nlocal][m]); } - m++; + m++; array[nlocal][0] = extra[nlocal][m++]; array[nlocal][1] = extra[nlocal][m++]; @@ -561,7 +561,7 @@ int FixSRP::size_restart(int nlocal) } /* ---------------------------------------------------------------------- - pack global state of Fix + pack global state of Fix ------------------------------------------------------------------------- */ void FixSRP::write_restart(FILE *fp) diff --git a/src/USER-MISC/fix_srp.h b/src/USER-MISC/fix_srp.h index e64629b6be..859b0020b0 100644 --- a/src/USER-MISC/fix_srp.h +++ b/src/USER-MISC/fix_srp.h @@ -43,7 +43,7 @@ class FixSRP : public Fix { int unpack_exchange(int, double *); int pack_border(int, int *, double *); int unpack_border(int, int, double *); - void post_run(); + void post_run(); int pack_restart(int, double*); void unpack_restart(int, int); diff --git a/src/USER-MISC/pair_srp.cpp b/src/USER-MISC/pair_srp.cpp index 0fcacb0790..fedbe0a302 100644 --- a/src/USER-MISC/pair_srp.cpp +++ b/src/USER-MISC/pair_srp.cpp @@ -10,7 +10,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ - + /* ---------------------------------------------------------------------- Contributing authors: Timothy Sirk (ARL), Pieter in't Veld (BASF) @@ -25,16 +25,16 @@ There is an example script for this package in examples/USER/srp. Please contact Timothy Sirk for questions (tim.sirk@us.army.mil). ------------------------------------------------------------------------- */ -#include "stdlib.h" -#include "pair_srp.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" -#include "domain.h" +#include "stdlib.h" +#include "pair_srp.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "domain.h" #include "modify.h" #include "fix.h" #include "fix_srp.h" @@ -42,9 +42,9 @@ Please contact Timothy Sirk for questions (tim.sirk@us.army.mil). #include "output.h" #include "string.h" #include "citeme.h" - -using namespace LAMMPS_NS; - + +using namespace LAMMPS_NS; + #define SMALL 1.0e-10 #define BIG 1e10 #define ONETWOBIT 0x40000000 @@ -65,10 +65,10 @@ static int srp_instance = 0; set size of pair comms in constructor ---------------------------------------------------------------------- */ -PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) +PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) { - writedata = 1; - + writedata = 1; + if (lmp->citeme) lmp->citeme->add(cite_srp); nextra = 1; @@ -77,8 +77,8 @@ PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) fix_id[0] = '0' + srp_instance / 10; fix_id[1] = '0' + srp_instance % 10; ++srp_instance; -} - +} + /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ @@ -102,19 +102,19 @@ void PairSRP::allocate() maxcount = 0; } -/* ---------------------------------------------------------------------- - free - ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + free + ------------------------------------------------------------------------- */ -PairSRP::~PairSRP() -{ - if (allocated) - { - memory->destroy(setflag); - memory->destroy(cutsq); - memory->destroy(cut); - memory->destroy(a0); - memory->destroy(segment); +PairSRP::~PairSRP() +{ + if (allocated) + { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cut); + memory->destroy(a0); + memory->destroy(segment); } // check nfix in case all fixes have already been deleted @@ -122,24 +122,24 @@ PairSRP::~PairSRP() free(fix_id); } -/* ---------------------------------------------------------------------- - compute bond-bond repulsions - ------------------------------------------------------------------------- */ - -void PairSRP::compute(int eflag, int vflag) +/* ---------------------------------------------------------------------- + compute bond-bond repulsions + ------------------------------------------------------------------------- */ + +void PairSRP::compute(int eflag, int vflag) { - // setup energy and virial - if (eflag || vflag) - ev_setup(eflag, vflag); - else - evflag = vflag_fdotr = 0; + // setup energy and virial + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; - double **x = atom->x; - double **f = atom->f; - int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - int i0, i1, j0, j1; + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + int i0, i1, j0, j1; int i,j,ii,jj,inum,jnum; double dijsq, dij; @@ -154,8 +154,8 @@ void PairSRP::compute(int eflag, int vflag) double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1; double fx, fy, fz; - // mapping global to local for atoms inside bond particles - // exclude 1-2 neighs if requested + // mapping global to local for atoms inside bond particles + // exclude 1-2 neighs if requested if (neighbor->ago == 0){ remapBonds(nall); if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh); @@ -163,9 +163,9 @@ void PairSRP::compute(int eflag, int vflag) // this pair style only used with hybrid // due to exclusions - // each atom i is type bptype - // each neigh j is type bptype - + // each atom i is type bptype + // each neigh j is type bptype + // using midpoint distance option if(midpoint){ @@ -193,18 +193,18 @@ void PairSRP::compute(int eflag, int vflag) // midpt dist bond 0 and 1 dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]); - dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]); - dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]); + dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]); + dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]); dijsq = dx*dx + dy*dy + dz*dz; if (dijsq < cutsq[bptype][bptype]){ dij = sqrt(dijsq); - if (dij < SMALL) + if (dij < SMALL) continue; // dij can be 0.0 with soft potentials wd = 1.0 - dij / cut[bptype][bptype]; - fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule + fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule // force for bond 0, beads 0,1 //force between bonds @@ -244,10 +244,10 @@ void PairSRP::compute(int eflag, int vflag) } } } - } + } else{ // using min distance option - + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -274,13 +274,13 @@ void PairSRP::compute(int eflag, int vflag) if (dijsq < cutsq[bptype][bptype]){ - dij = sqrt(dijsq); + dij = sqrt(dijsq); if (dij < SMALL) continue; // dij can be 0.0 with soft potentials wd = 1.0 - dij / cut[bptype][bptype]; - fpair = a0[bptype][bptype] * wd / dij; + fpair = a0[bptype][bptype] * wd / dij; // force for bond 0, beads 0,1 lever0 = 0.5 + ti; // assign force according to lever rule @@ -364,7 +364,7 @@ void PairSRP::settings(int narg, char **arg) int iarg = 3; // default exclude 1-2 // scaling for 1-2, etc not supported - exclude = 1; + exclude = 1; while (iarg < narg) { if (strcmp(arg[iarg],"exclude") == 0) { @@ -392,7 +392,7 @@ void PairSRP::settings(int narg, char **arg) } /* ---------------------------------------------------------------------- - set coeffs + set coeffs ------------------------------------------------------------------------- */ void PairSRP::coeff(int narg, char **arg) @@ -459,7 +459,7 @@ void PairSRP::init_style() // bond particles do not contribute to energy or virial // bond particles do not belong to group all - // but thermo normalization is by nall + // but thermo normalization is by nall // therefore should turn off normalization int me; MPI_Comm_rank(world,&me); @@ -467,7 +467,7 @@ void PairSRP::init_style() arg1[0] = (char *) "norm"; arg1[1] = (char *) "no"; output->thermo->modify_params(2, arg1); - if (me == 0) + if (me == 0) error->message(FLERR,"Thermo normalization turned off by pair srp"); neighbor->request(this,instance_me); @@ -487,13 +487,13 @@ double PairSRP::init_one(int i, int j) return cut[i][j]; } - -/* ---------------------------------------------------------------------- + +/* ---------------------------------------------------------------------- find min distance for bonds i0/j0 and i1/j1 - ------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1) -{ - // move these outside the loop +{ + // move these outside the loop double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj; double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1; @@ -503,7 +503,7 @@ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, diffz0 = x[j0][2] - x[i0][2]; // compute midpt dist from 1st atom, 2nd bond - diffx1 = x[j1][0] - x[i1][0]; + diffx1 = x[j1][0] - x[i1][0]; diffy1 = x[j1][1] - x[i1][1]; diffz1 = x[j1][2] - x[i1][2]; @@ -519,22 +519,22 @@ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, denom = RiRj*RiRj - RiRi*RjRj; // handle case of parallel lines - // reduce to midpt distance - if (fabs(denom) < SMALL){ + // reduce to midpt distance + if (fabs(denom) < SMALL){ if(denom < 0) denom = -BIG; else denom = BIG; - } + } - // calc ti + // calc ti termx0 = RiRj*diffx1 - RjRj*diffx0; termy0 = RiRj*diffy1 - RjRj*diffy0; termz0 = RiRj*diffz1 - RjRj*diffz0; num0 = dPx*termx0 + dPy*termy0 + dPz*termz0; ti = num0 / denom; - if (ti > 0.5) ti = 0.5; - if (ti < -0.5) ti = -0.5; - - // calc tj + if (ti > 0.5) ti = 0.5; + if (ti < -0.5) ti = -0.5; + + // calc tj termx1 = RiRj*diffx0 - RiRi*diffx1; termy1 = RiRj*diffy0 - RiRi*diffy1; termz1 = RiRj*diffz0 - RiRi*diffz1; @@ -542,28 +542,28 @@ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, tj = -num1/ denom; if (tj > 0.5) tj = 0.5; if (tj < -0.5) tj = -0.5; - - // min dist + + // min dist dx = dPx - ti*diffx0 + tj*diffx1; dy = dPy - ti*diffy0 + tj*diffy1; dz = dPz - ti*diffz0 + tj*diffz1; -} +} -/* -------------------------------------------------------- +/* -------------------------------------------------------- map global id of atoms in stored by each bond particle - ------------------------------------------------------- */ + ------------------------------------------------------- */ inline void PairSRP::remapBonds(int &nall) { if(nall > maxcount){ - memory->grow(segment, nall, 2, "pair:segment"); + memory->grow(segment, nall, 2, "pair:segment"); maxcount = nall; } // loop over all bond particles // each bond paricle holds two bond atoms // map global ids of bond atoms to local ids - // might not be able to map both bond atoms of j, if j is outside neighcut - // these are not on neighlist, so are not used + // might not be able to map both bond atoms of j, if j is outside neighcut + // these are not on neighlist, so are not used int tmp; srp = f_srp->array_atom; @@ -580,18 +580,18 @@ inline void PairSRP::remapBonds(int &nall) } } -/* -------------------------------------------------------- +/* -------------------------------------------------------- add exclusions for 1-2 neighs, if requested -more complex exclusions or scaling probably not needed - ------------------------------------------------------- */ +more complex exclusions or scaling probably not needed + ------------------------------------------------------- */ inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh) { int i0, i1, j0, j1; int i,j,ii,jj,jnum; // encode neighs with exclusions - // only need 1-2 info for normal uses of srp - // add 1-3, etc later if ever needed + // only need 1-2 info for normal uses of srp + // add 1-3, etc later if ever needed for (ii = 0; ii < inum; ii++) { @@ -610,7 +610,7 @@ inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &nu i1 = segment[j][0]; j1 = segment[j][1]; - // check for a 1-2 neigh + // check for a 1-2 neigh if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){ j |= ONETWOBIT; jlist[jj] = j; @@ -639,7 +639,7 @@ void PairSRP::write_data_all(FILE *fp) for (int j = i; j <= atom->ntypes; j++) fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]); } - + /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ diff --git a/src/USER-MISC/pair_srp.h b/src/USER-MISC/pair_srp.h index 95cc5a5c94..ecd352cdeb 100644 --- a/src/USER-MISC/pair_srp.h +++ b/src/USER-MISC/pair_srp.h @@ -44,9 +44,9 @@ class PairSRP : public Pair { inline void onetwoexclude(int* &, int &, int* &, int* &, int** &); inline void remapBonds(int &); void allocate(); - void getMinDist(double** &, double &, double &, double &, double &, + void getMinDist(double** &, double &, double &, double &, double &, double &, int &, int &, int &, int &); - bool min, midpoint; + bool min, midpoint; double **cut; double **a0; double **srp; diff --git a/src/USER-OMP/fix_nh_asphere_omp.cpp b/src/USER-OMP/fix_nh_asphere_omp.cpp index 4c23c7b14d..c9993e50e1 100644 --- a/src/USER-OMP/fix_nh_asphere_omp.cpp +++ b/src/USER-OMP/fix_nh_asphere_omp.cpp @@ -79,7 +79,7 @@ void FixNHAsphereOMP::nve_v() const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal; int i; - // standard nve_v velocity update. for efficiency the loop is + // standard nve_v velocity update. for efficiency the loop is // merged with FixNHOMP instead of calling it for the COM update. #if defined(_OPENMP) @@ -122,7 +122,7 @@ void FixNHAsphereOMP::nve_x() // update quaternion a full step via Richardson iteration // returns new normalized quaternion // principal moments of inertia - + #if defined(_OPENMP) #pragma omp parallel for default(none) private(i) schedule(static) #endif diff --git a/src/USER-OMP/fix_nh_sphere_omp.cpp b/src/USER-OMP/fix_nh_sphere_omp.cpp index 739124575e..98a0141bd1 100644 --- a/src/USER-OMP/fix_nh_sphere_omp.cpp +++ b/src/USER-OMP/fix_nh_sphere_omp.cpp @@ -81,7 +81,7 @@ void FixNHSphereOMP::nve_v() const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal; int i; - // standard nve_v velocity update. for efficiency the loop is + // standard nve_v velocity update. for efficiency the loop is // merged with FixNHOMP instead of calling it for the COM update. // update omega for all particles diff --git a/src/USER-OMP/fix_rigid_nh_omp.cpp b/src/USER-OMP/fix_rigid_nh_omp.cpp index a35d849038..0e347eb128 100644 --- a/src/USER-OMP/fix_rigid_nh_omp.cpp +++ b/src/USER-OMP/fix_rigid_nh_omp.cpp @@ -70,7 +70,7 @@ void FixRigidNHOMP::initial_integrate(int vflag) scale_t[0] = scale_t[1] = scale_t[2] = tmp; tmp = exp(-dtq * eta_dot_r[0]); scale_r = tmp; - } + } if (pstat_flag) { akin_t = akin_r = 0.0; @@ -86,7 +86,7 @@ void FixRigidNHOMP::initial_integrate(int vflag) tmp = dtq * epsilon_dot[2]; scale_v[2] = dtv * exp(tmp) * maclaurin_series(tmp); } - + // update xcm, vcm, quat, conjqm and angmom double akt=0.0, akr=0.0; int ibody; @@ -97,24 +97,24 @@ void FixRigidNHOMP::initial_integrate(int vflag) for (ibody = 0; ibody < nbody; ibody++) { double mbody[3],tbody[3],fquat[4]; const double dtf2 = dtf * 2.0; - + // step 1.1 - update vcm by 1/2 step - + const double dtfm = dtf / masstotal[ibody]; vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; - + if (tstat_flag || pstat_flag) { vcm[ibody][0] *= scale_t[0]; vcm[ibody][1] *= scale_t[1]; vcm[ibody][2] *= scale_t[2]; - + double tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + vcm[ibody][2]*vcm[ibody][2]; akt += masstotal[ibody]*tmp; } - + // step 1.2 - update xcm by full step if (!pstat_flag) { @@ -126,29 +126,29 @@ void FixRigidNHOMP::initial_integrate(int vflag) xcm[ibody][1] += scale_v[1] * vcm[ibody][1]; xcm[ibody][2] += scale_v[2] * vcm[ibody][2]; } - + // step 1.3 - apply torque (body coords) to quaternion momentum - + torque[ibody][0] *= tflag[ibody][0]; torque[ibody][1] *= tflag[ibody][1]; torque[ibody][2] *= tflag[ibody][2]; - + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], torque[ibody],tbody); MathExtra::quatvec(quat[ibody],tbody,fquat); - + conjqm[ibody][0] += dtf2 * fquat[0]; conjqm[ibody][1] += dtf2 * fquat[1]; conjqm[ibody][2] += dtf2 * fquat[2]; conjqm[ibody][3] += dtf2 * fquat[3]; - + if (tstat_flag || pstat_flag) { conjqm[ibody][0] *= scale_r; conjqm[ibody][1] *= scale_r; conjqm[ibody][2] *= scale_r; conjqm[ibody][3] *= scale_r; } - + // step 1.4 to 1.13 - use no_squish rotate to update p and q MathExtra::no_squish_rotate(3,conjqm[ibody],quat[ibody],inertia[ibody],dtq); @@ -160,22 +160,22 @@ void FixRigidNHOMP::initial_integrate(int vflag) // update exyz_space // transform p back to angmom // update angular velocity - + MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody]); MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], mbody,angmom[ibody]); - + angmom[ibody][0] *= 0.5; angmom[ibody][1] *= 0.5; angmom[ibody][2] *= 0.5; - + MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); - + if (tstat_flag || pstat_flag) { - akr += angmom[ibody][0]*omega[ibody][0] + + akr += angmom[ibody][0]*omega[ibody][0] + angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; } } // end of parallel for @@ -205,14 +205,14 @@ void FixRigidNHOMP::initial_integrate(int vflag) if (vflag) v_setup(vflag); else evflag = 0; - + // remap simulation box by 1/2 step if (pstat_flag) remap(); - + // set coords/orient and velocity/rotation of atoms in rigid bodies // from quarternion and omega - + if (triclinic) if (evflag) set_xv_thr<1,1>(); @@ -230,7 +230,7 @@ void FixRigidNHOMP::initial_integrate(int vflag) if (pstat_flag) { remap(); if (kspace_flag) force->kspace->setup(); - } + } } /* ---------------------------------------------------------------------- */ @@ -240,7 +240,7 @@ void FixRigidNHOMP::final_integrate() double scale_t[3],scale_r; // compute scale variables - + scale_t[0] = scale_t[1] = scale_t[2] = 1.0; scale_r = 1.0; @@ -248,17 +248,17 @@ void FixRigidNHOMP::final_integrate() double tmp = exp(-1.0 * dtq * eta_dot_t[0]); scale_t[0] = scale_t[1] = scale_t[2] = tmp; scale_r = exp(-1.0 * dtq * eta_dot_r[0]); - } - + } + if (pstat_flag) { scale_t[0] *= exp(-dtq * (epsilon_dot[0] + mtk_term2)); scale_t[1] *= exp(-dtq * (epsilon_dot[1] + mtk_term2)); scale_t[2] *= exp(-dtq * (epsilon_dot[2] + mtk_term2)); scale_r *= exp(-dtq * (pdim * mtk_term2)); - + akin_t = akin_r = 0.0; } - + double * const * _noalias const x = atom->x; const dbl3_t * _noalias const f = (dbl3_t *) atom->f[0]; const double * const * const torque_one = atom->torque; @@ -306,7 +306,7 @@ void FixRigidNHOMP::final_integrate() // we likely have only a rather number of groups so we loops // over bodies and thread over all atoms for each of them. - + for (int ib = 0; ib < nbody; ++ib) { double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; int i; @@ -361,7 +361,7 @@ void FixRigidNHOMP::final_integrate() #else const int tid = 0; #endif - + for (int i = 0; i < nlocal; i++) { const int ibody = body[i]; if ((ibody < 0) || (ibody % nthreads != tid)) continue; @@ -423,28 +423,28 @@ void FixRigidNHOMP::final_integrate() vcm[ibody][1] *= scale_t[1]; vcm[ibody][2] *= scale_t[2]; } - + vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0]; vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1]; vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2]; - + if (pstat_flag) { double tmp = vcm[ibody][0]*vcm[ibody][0] + vcm[ibody][1]*vcm[ibody][1] + vcm[ibody][2]*vcm[ibody][2]; akt += masstotal[ibody]*tmp; } - + // update conjqm, then transform to angmom, set velocity again // virial is already setup from initial_integrate - + torque[ibody][0] *= tflag[ibody][0]; torque[ibody][1] *= tflag[ibody][1]; torque[ibody][2] *= tflag[ibody][2]; - + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody], ez_space[ibody],torque[ibody],tbody); MathExtra::quatvec(quat[ibody],tbody,fquat); - + if (tstat_flag || pstat_flag) { conjqm[ibody][0] = scale_r * conjqm[ibody][0] + dtf2 * fquat[0]; conjqm[ibody][1] = scale_r * conjqm[ibody][1] + dtf2 * fquat[1]; @@ -460,17 +460,17 @@ void FixRigidNHOMP::final_integrate() MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], mbody,angmom[ibody]); - + angmom[ibody][0] *= 0.5; angmom[ibody][1] *= 0.5; - angmom[ibody][2] *= 0.5; - + angmom[ibody][2] *= 0.5; + MathExtra::angmom_to_omega(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); - + if (pstat_flag) { - akr += angmom[ibody][0]*omega[ibody][0] + - angmom[ibody][1]*omega[ibody][1] + + akr += angmom[ibody][0]*omega[ibody][0] + + angmom[ibody][1]*omega[ibody][1] + angmom[ibody][2]*omega[ibody][2]; } } // end of parallel for @@ -523,11 +523,11 @@ void FixRigidNHOMP::remap() const int nlocal = atom->nlocal; // epsilon is not used, except for book-keeping - + for (int i = 0; i < 3; i++) epsilon[i] += dtq * epsilon_dot[i]; - + // convert pertinent atoms and rigid bodies to lamda coords - + if (allremap) domain->x2lamda(nlocal); else { int i; @@ -538,13 +538,13 @@ void FixRigidNHOMP::remap() if (mask[i] & dilate_group_bit) domain->x2lamda(x[i],x[i]); } - + if (nrigid) for (int i = 0; i < nrigidfix; i++) modify->fix[rfix[i]]->deform(0); - + // reset global and local box to new size/shape - + for (int i = 0; i < 3; i++) { if (p_flag[i]) { const double oldlo = domain->boxlo[i]; @@ -558,9 +558,9 @@ void FixRigidNHOMP::remap() domain->set_global_box(); domain->set_local_box(); - + // convert pertinent atoms and rigid bodies back to box coords - + if (allremap) domain->lamda2x(nlocal); else { int i; @@ -901,7 +901,7 @@ void FixRigidNHOMP::set_v_thr() } } - // set omega, angmom of each extended particle + // set omega, angmom of each extended particle // XXX: extended particle info not yet multi-threaded if (extended) { diff --git a/src/USER-OMP/fix_rigid_nh_omp.h b/src/USER-OMP/fix_rigid_nh_omp.h index 9210d7e5dc..f2902ca333 100644 --- a/src/USER-OMP/fix_rigid_nh_omp.h +++ b/src/USER-OMP/fix_rigid_nh_omp.h @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -26,7 +26,7 @@ class FixRigidNHOMP : public FixRigidNH { virtual void initial_integrate(int); virtual void final_integrate(); - virtual void remap(); + virtual void remap(); private: // copied from FixRigidOMP template void set_xv_thr(); diff --git a/src/USER-OMP/fix_rigid_nph_omp.cpp b/src/USER-OMP/fix_rigid_nph_omp.cpp index 99506bbebe..bd2ba5f6fb 100644 --- a/src/USER-OMP/fix_rigid_nph_omp.cpp +++ b/src/USER-OMP/fix_rigid_nph_omp.cpp @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -29,24 +29,24 @@ using namespace LAMMPS_NS; FixRigidNPHOMP::FixRigidNPHOMP(LAMMPS *lmp, int narg, char **arg) : FixRigidNHOMP(lmp, narg, arg) -{ +{ // other setting are made by parent scalar_flag = 1; restart_global = 1; box_change_size = 1; extscalar = 1; - + // error checks if (pstat_flag == 0) error->all(FLERR,"Pressure control must be used with fix nph/omp"); if (tstat_flag == 1) error->all(FLERR,"Temperature control must not be used with fix nph/omp"); - if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || + if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) error->all(FLERR,"Target pressure for fix rigid/nph/omp cannot be 0.0"); - + // convert input periods to frequency p_freq[0] = p_freq[1] = p_freq[2] = 0.0; @@ -67,15 +67,15 @@ FixRigidNPHOMP::FixRigidNPHOMP(LAMMPS *lmp, int narg, char **arg) : char **newarg = new char*[3]; newarg[0] = id_temp; newarg[1] = (char *) "all"; - newarg[2] = (char *) "temp"; + newarg[2] = (char *) "temp"; modify->add_compute(3,newarg); delete [] newarg; tcomputeflag = 1; - + // create a new compute pressure style // id = fix-ID + press, compute group = all // pass id_temp as 4th arg to pressure constructor - + n = strlen(id) + 7; id_press = new char[n]; strcpy(id_press,id); diff --git a/src/USER-OMP/fix_rigid_nph_omp.h b/src/USER-OMP/fix_rigid_nph_omp.h index f64a743fea..31b95fa02f 100644 --- a/src/USER-OMP/fix_rigid_nph_omp.h +++ b/src/USER-OMP/fix_rigid_nph_omp.h @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. diff --git a/src/USER-OMP/fix_rigid_npt_omp.cpp b/src/USER-OMP/fix_rigid_npt_omp.cpp index f68b68c29f..5d41d5f296 100644 --- a/src/USER-OMP/fix_rigid_npt_omp.cpp +++ b/src/USER-OMP/fix_rigid_npt_omp.cpp @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -29,21 +29,21 @@ using namespace LAMMPS_NS; FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) : FixRigidNHOMP(lmp, narg, arg) -{ +{ // other setting are made by parent scalar_flag = 1; restart_global = 1; box_change_size = 1; extscalar = 1; - + // error checks if (tstat_flag == 0 || pstat_flag == 0) error->all(FLERR,"Did not set temp or press for fix rigid/npt/omp"); if (t_start <= 0.0 || t_stop <= 0.0) error->all(FLERR,"Target temperature for fix rigid/npt/omp cannot be 0.0"); - if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || + if (p_start[0] < 0.0 || p_start[1] < 0.0 || p_start[2] < 0.0 || p_stop[0] < 0.0 || p_stop[1] < 0.0 || p_stop[2] < 0.0) error->all(FLERR,"Target pressure for fix rigid/npt/omp cannot be 0.0"); @@ -53,7 +53,7 @@ FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) : if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); - if (t_order != 3 && t_order != 5) + if (t_order != 3 && t_order != 5) error->all(FLERR,"Fix_modify order must be 3 or 5"); // convert input periods to frequency @@ -79,15 +79,15 @@ FixRigidNPTOMP::FixRigidNPTOMP(LAMMPS *lmp, int narg, char **arg) : char **newarg = new char*[3]; newarg[0] = id_temp; newarg[1] = (char *) "all"; - newarg[2] = (char *) "temp"; + newarg[2] = (char *) "temp"; modify->add_compute(3,newarg); delete [] newarg; tcomputeflag = 1; - + // create a new compute pressure style // id = fix-ID + press, compute group = all // pass id_temp as 4th arg to pressure constructor - + n = strlen(id) + 7; id_press = new char[n]; strcpy(id_press,id); diff --git a/src/USER-OMP/fix_rigid_npt_omp.h b/src/USER-OMP/fix_rigid_npt_omp.h index 627b5a8c29..4f17d0ba55 100644 --- a/src/USER-OMP/fix_rigid_npt_omp.h +++ b/src/USER-OMP/fix_rigid_npt_omp.h @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. diff --git a/src/USER-OMP/fix_rigid_nve_omp.cpp b/src/USER-OMP/fix_rigid_nve_omp.cpp index 51ed34fae8..1da3c29fa7 100644 --- a/src/USER-OMP/fix_rigid_nve_omp.cpp +++ b/src/USER-OMP/fix_rigid_nve_omp.cpp @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. diff --git a/src/USER-OMP/fix_rigid_nve_omp.h b/src/USER-OMP/fix_rigid_nve_omp.h index eea45ee358..464dde528b 100644 --- a/src/USER-OMP/fix_rigid_nve_omp.h +++ b/src/USER-OMP/fix_rigid_nve_omp.h @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. diff --git a/src/USER-OMP/fix_rigid_nvt_omp.cpp b/src/USER-OMP/fix_rigid_nvt_omp.cpp index 758dcce9bf..5e367899d9 100644 --- a/src/USER-OMP/fix_rigid_nvt_omp.cpp +++ b/src/USER-OMP/fix_rigid_nvt_omp.cpp @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. @@ -26,13 +26,13 @@ using namespace LAMMPS_NS; FixRigidNVTOMP::FixRigidNVTOMP(LAMMPS *lmp, int narg, char **arg) : FixRigidNHOMP(lmp, narg, arg) -{ +{ // other settings are made by parent scalar_flag = 1; restart_global = 1; extscalar = 1; - + // error checking // convert input period to frequency @@ -46,6 +46,6 @@ FixRigidNVTOMP::FixRigidNVTOMP(LAMMPS *lmp, int narg, char **arg) : if (t_chain < 1) error->all(FLERR,"Illegal fix_modify command"); if (t_iter < 1) error->all(FLERR,"Illegal fix_modify command"); - if (t_order != 3 && t_order != 5) - error->all(FLERR,"Fix_modify order must be 3 or 5"); + if (t_order != 3 && t_order != 5) + error->all(FLERR,"Fix_modify order must be 3 or 5"); } diff --git a/src/USER-OMP/fix_rigid_nvt_omp.h b/src/USER-OMP/fix_rigid_nvt_omp.h index e1583aa4df..d083834d82 100644 --- a/src/USER-OMP/fix_rigid_nvt_omp.h +++ b/src/USER-OMP/fix_rigid_nvt_omp.h @@ -5,7 +5,7 @@ Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under + certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. diff --git a/src/USER-OMP/fix_rigid_omp.cpp b/src/USER-OMP/fix_rigid_omp.cpp index 145df99f31..8e069abfbf 100644 --- a/src/USER-OMP/fix_rigid_omp.cpp +++ b/src/USER-OMP/fix_rigid_omp.cpp @@ -155,9 +155,9 @@ void FixRigidOMP::final_integrate() } else if (rstyle == GROUP) { - // we likely have only a rather number of groups so we loop + // we likely have only a rather number of groups so we loops // over bodies and thread over all atoms for each of them. - + for (int ib = 0; ib < nbody; ++ib) { double s0=0.0,s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0; int i; @@ -212,7 +212,7 @@ void FixRigidOMP::final_integrate() #else const int tid = 0; #endif - + for (int i = 0; i < nlocal; i++) { const int ibody = body[i]; if ((ibody < 0) || (ibody % nthreads != tid)) continue; @@ -616,7 +616,7 @@ void FixRigidOMP::set_v_thr() } } - // set omega, angmom of each extended particle + // set omega, angmom of each extended particle // XXX: extended particle info not yet multi-threaded if (extended) { diff --git a/src/USER-OMP/fix_rigid_small_omp.cpp b/src/USER-OMP/fix_rigid_small_omp.cpp index a260899aef..e419bc5227 100644 --- a/src/USER-OMP/fix_rigid_small_omp.cpp +++ b/src/USER-OMP/fix_rigid_small_omp.cpp @@ -120,7 +120,7 @@ void FixRigidSmallOMP::final_integrate() const int nlocal = atom->nlocal; const int nthreads=comm->nthreads; int i, ibody; - + #if defined(_OPENMP) #pragma omp parallel for default(none) private(ibody) schedule(static) #endif @@ -139,7 +139,7 @@ void FixRigidSmallOMP::final_integrate() #if defined(_OPENMP) #pragma omp parallel default(none) private(i,ibody) #endif - { + { #if defined(_OPENMP) const int tid = omp_get_thread_num(); #else @@ -149,7 +149,7 @@ void FixRigidSmallOMP::final_integrate() for (i = 0; i < nlocal; i++) { ibody = atom2body[i]; if ((ibody < 0) || (ibody % nthreads != tid)) continue; - + Body &b = body[ibody]; double unwrap[3]; @@ -203,7 +203,7 @@ void FixRigidSmallOMP::final_integrate() } // update vcm and angmom, recompute omega - + #if defined(_OPENMP) #pragma omp parallel for default(none) private(ibody) schedule(static) #endif @@ -562,7 +562,7 @@ void FixRigidSmallOMP::set_v_thr() } } - // set omega, angmom of each extended particle + // set omega, angmom of each extended particle // XXX: extended particle info not yet multi-threaded if (extended) { diff --git a/src/USER-OMP/thr_omp.cpp b/src/USER-OMP/thr_omp.cpp index e636289848..713e7fdf5d 100644 --- a/src/USER-OMP/thr_omp.cpp +++ b/src/USER-OMP/thr_omp.cpp @@ -22,7 +22,6 @@ #include "force.h" #include "memory.h" #include "modify.h" -#include "compute.h" #include "neighbor.h" #include "timer.h" diff --git a/src/USER-PHONON/fix_phonon.cpp b/src/USER-PHONON/fix_phonon.cpp index 9d601d17bf..b0d09c146f 100644 --- a/src/USER-PHONON/fix_phonon.cpp +++ b/src/USER-PHONON/fix_phonon.cpp @@ -666,10 +666,9 @@ void FixPhonon::postprocess( ) } // to get Phi = KT.G^-1; normalization of FFTW data is done here - double boltz = force->boltz, *kbtsqrt, TempAve = 0.; + double boltz = force->boltz, kbtsqrt[sysdim], TempAve = 0.; double TempFac = inv_neval*inv_nTemp; double NormFac = TempFac*double(ntotal); - kbtsqrt = new double[sysdim]; for (idim = 0; idim < sysdim; ++idim){ kbtsqrt[idim] = sqrt(TempSum[idim]*NormFac); @@ -683,7 +682,6 @@ void FixPhonon::postprocess( ) for (idim = 0; idim < fft_dim; ++idim) for (jdim = 0; jdim < fft_dim; ++jdim) Phi_q[idq][ndim++] *= kbtsqrt[idim%sysdim]*kbtsqrt[jdim%sysdim]; } - delete[] kbtsqrt; // to collect all local Phi_q to root displs[0]=0; @@ -692,8 +690,7 @@ void FixPhonon::postprocess( ) MPI_Gatherv(Phi_q[0],mynq*fft_dim2*2,MPI_DOUBLE,Phi_all[0],recvcnts,displs,MPI_DOUBLE,0,world); // to collect all basis info and averaged it on root - double *basis_root; - basis_root = new double[fft_dim]; + double basis_root[fft_dim]; if (fft_dim > sysdim) MPI_Reduce(&basis[1][0], &basis_root[sysdim], fft_dim-sysdim, MPI_DOUBLE, MPI_SUM, 0, world); if (me == 0){ // output dynamic matrix by root @@ -733,7 +730,6 @@ void FixPhonon::postprocess( ) fwrite(M_inv_sqrt, sizeof(double),nucell, fp_bin); fclose(fp_bin); - delete[] basis_root; // write log file, here however, it is the dynamical matrix that is written for (int i = 0; i < 60; ++i) fprintf(flog,"#"); fprintf(flog,"\n"); diff --git a/src/atom.cpp b/src/atom.cpp index fc7641f9ba..23f43c2e00 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -735,23 +735,23 @@ void Atom::deallocate_topology() memory->destroy(atom->angle_atom3); atom->angle_type = NULL; atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = NULL; - + memory->destroy(atom->dihedral_type); memory->destroy(atom->dihedral_atom1); memory->destroy(atom->dihedral_atom2); memory->destroy(atom->dihedral_atom3); memory->destroy(atom->dihedral_atom4); atom->dihedral_type = NULL; - atom->dihedral_atom1 = atom->dihedral_atom2 = + atom->dihedral_atom1 = atom->dihedral_atom2 = atom->dihedral_atom3 = atom->dihedral_atom4 = NULL; - + memory->destroy(atom->improper_type); memory->destroy(atom->improper_atom1); memory->destroy(atom->improper_atom2); memory->destroy(atom->improper_atom3); memory->destroy(atom->improper_atom4); atom->improper_type = NULL; - atom->improper_atom1 = atom->improper_atom2 = + atom->improper_atom1 = atom->improper_atom2 = atom->improper_atom3 = atom->improper_atom4 = NULL; } @@ -865,7 +865,7 @@ void Atom::data_atoms(int n, char *buf, tagint id_offset, int type_offset, (((imageint) (atoi(values[iptr+2]) + IMGMAX) & IMGMASK) << IMG2BITS); else imagedata = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; - + xdata[0] = atof(values[xptr]); xdata[1] = atof(values[xptr+1]); xdata[2] = atof(values[xptr+2]); @@ -1078,7 +1078,7 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset) for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; - sscanf(buf,"%d %d " + sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); if (id_offset) { @@ -1163,7 +1163,7 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset) for (int i = 0; i < n; i++) { next = strchr(buf,'\n'); *next = '\0'; - sscanf(buf,"%d %d " + sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, &tmp,&itype,&atom1,&atom2,&atom3,&atom4); if (id_offset) { @@ -1540,7 +1540,7 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom, if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom]; if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom]; if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom]; - else if (rmass_flag) + else if (rmass_flag) rmass[ilocal] = 4.0*MY_PI/3.0 * radius[ilocal]*radius[ilocal]*radius[ilocal]; diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index c348bc806f..b4b1d3f210 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -65,8 +65,8 @@ AtomVecBody::~AtomVecBody() { int nall = nlocal_bonus + nghost_bonus; for (int i = 0; i < nall; i++) { - bptr->icp->put(bonus[i].iindex); - bptr->dcp->put(bonus[i].dindex); + icp->put(bonus[i].iindex); + dcp->put(bonus[i].dindex); } memory->sfree(bonus); @@ -95,6 +95,8 @@ void AtomVecBody::process_args(int narg, char **arg) else error->all(FLERR,"Unknown body style"); bptr->avec = this; + icp = bptr->icp; + dcp = bptr->dcp; // max size of forward/border comm // 7,16 are packed in pack_comm/pack_border @@ -189,8 +191,8 @@ void AtomVecBody::copy(int i, int j, int delflag) // if deleting atom J via delflag and J has bonus data, then delete it if (delflag && body[j] >= 0) { - bptr->icp->put(bonus[body[j]].iindex); - bptr->dcp->put(bonus[body[j]].dindex); + icp->put(bonus[body[j]].iindex); + dcp->put(bonus[body[j]].dindex); copy_bonus(nlocal_bonus-1,body[j]); nlocal_bonus--; } @@ -226,8 +228,8 @@ void AtomVecBody::clear_bonus() { int nall = nlocal_bonus + nghost_bonus; for (int i = nlocal_bonus; i < nall; i++) { - bptr->icp->put(bonus[i].iindex); - bptr->dcp->put(bonus[i].dindex); + icp->put(bonus[i].iindex); + dcp->put(bonus[i].dindex); } nghost_bonus = 0; } @@ -828,8 +830,8 @@ void AtomVecBody::unpack_border(int n, int first, double *buf) inertia[2] = buf[m++]; bonus[j].ninteger = (int) ubuf(buf[m++]).i; bonus[j].ndouble = (int) ubuf(buf[m++]).i; - bonus[j].ivalue = bptr->icp->get(bonus[j].ninteger,bonus[j].iindex); - bonus[j].dvalue = bptr->dcp->get(bonus[j].ndouble,bonus[j].dindex); + bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex); + bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex); m += bptr->unpack_border_body(&bonus[j],&buf[m]); bonus[j].ilocal = i; body[i] = j; @@ -876,8 +878,8 @@ void AtomVecBody::unpack_border_vel(int n, int first, double *buf) inertia[2] = buf[m++]; bonus[j].ninteger = (int) ubuf(buf[m++]).i; bonus[j].ndouble = (int) ubuf(buf[m++]).i; - bonus[j].ivalue = bptr->icp->get(bonus[j].ninteger,bonus[j].iindex); - bonus[j].dvalue = bptr->dcp->get(bonus[j].ndouble,bonus[j].dindex); + bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex); + bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex); m += bptr->unpack_border_body(&bonus[j],&buf[m]); bonus[j].ilocal = i; body[i] = j; @@ -923,8 +925,8 @@ int AtomVecBody::unpack_border_hybrid(int n, int first, double *buf) inertia[2] = buf[m++]; bonus[j].ninteger = (int) ubuf(buf[m++]).i; bonus[j].ndouble = (int) ubuf(buf[m++]).i; - bonus[j].ivalue = bptr->icp->get(bonus[j].ninteger,bonus[j].iindex); - bonus[j].dvalue = bptr->dcp->get(bonus[j].ndouble,bonus[j].dindex); + bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex); + bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex); m += bptr->unpack_border_body(&bonus[j],&buf[m]); bonus[j].ilocal = i; body[i] = j; @@ -1027,9 +1029,9 @@ int AtomVecBody::unpack_exchange(double *buf) inertia[2] = buf[m++]; bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i; bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i; - bonus[nlocal_bonus].ivalue = bptr->icp->get(bonus[nlocal_bonus].ninteger, + bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger, bonus[nlocal_bonus].iindex); - bonus[nlocal_bonus].dvalue = bptr->dcp->get(bonus[nlocal_bonus].ndouble, + bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble, bonus[nlocal_bonus].dindex); memcpy(bonus[nlocal_bonus].ivalue,&buf[m], bonus[nlocal_bonus].ninteger*sizeof(int)); @@ -1179,9 +1181,9 @@ int AtomVecBody::unpack_restart(double *buf) inertia[2] = buf[m++]; bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i; bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i; - bonus[nlocal_bonus].ivalue = bptr->icp->get(bonus[nlocal_bonus].ninteger, + bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger, bonus[nlocal_bonus].iindex); - bonus[nlocal_bonus].dvalue = bptr->dcp->get(bonus[nlocal_bonus].ndouble, + bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble, bonus[nlocal_bonus].dindex); memcpy(bonus[nlocal_bonus].ivalue,&buf[m], bonus[nlocal_bonus].ninteger*sizeof(int)); @@ -1469,7 +1471,7 @@ bigint AtomVecBody::memory_usage() if (atom->memcheck("body")) bytes += memory->usage(body,nmax); bytes += nmax_bonus*sizeof(Bonus); - bytes += bptr->icp->size + bptr->dcp->size; + bytes += icp->size + dcp->size; int nall = nlocal_bonus + nghost_bonus; for (int i = 0; i < nall; i++) { diff --git a/src/delete_bonds.cpp b/src/delete_bonds.cpp index 6ed4d211d1..c4823cc01b 100644 --- a/src/delete_bonds.cpp +++ b/src/delete_bonds.cpp @@ -25,8 +25,6 @@ #include "special.h" #include "error.h" -#include - using namespace LAMMPS_NS; enum{MULTI,ATOM,BOND,ANGLE,DIHEDRAL,IMPROPER,STATS}; diff --git a/src/error.h b/src/error.h index 09ea0b0447..9e8a6491c2 100644 --- a/src/error.h +++ b/src/error.h @@ -21,7 +21,6 @@ namespace LAMMPS_NS { class Error : protected Pointers { public: Error(class LAMMPS *); - virtual ~Error() {} void universe_all(const char *, int, const char *); void universe_one(const char *, int, const char *); diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index 3fd807e373..521f0956bd 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -413,7 +413,7 @@ FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) : maxvar = 0; varatom = NULL; - maxatom = -1; + maxatom = 0; bin = NULL; nbins = maxbin = 0; diff --git a/src/fix_heat.cpp b/src/fix_heat.cpp index 746801f846..de23927fe9 100644 --- a/src/fix_heat.cpp +++ b/src/fix_heat.cpp @@ -81,7 +81,7 @@ FixHeat::FixHeat(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) scale = 1.0; - maxatom = -1; + maxatom = 0; vheat = NULL; vscale = NULL; } diff --git a/src/fix_nh.h b/src/fix_nh.h index 0e676b4a0d..71ea86eebc 100644 --- a/src/fix_nh.h +++ b/src/fix_nh.h @@ -120,11 +120,10 @@ class FixNH : public Fix { double fixedpoint[3]; // location of dilation fixed-point void couple(); + void remap(); void nhc_temp_integrate(); void nhc_press_integrate(); - virtual void remap(); - virtual void nve_x(); // may be overwritten by child classes virtual void nve_v(); virtual void nh_v_press(); diff --git a/src/force.h b/src/force.h index 81371ee451..a3fca51020 100644 --- a/src/force.h +++ b/src/force.h @@ -75,7 +75,7 @@ class Force : protected Pointers { int special_extra; // extra space for added bonds Force(class LAMMPS *); - virtual ~Force(); + ~Force(); void init(); void create_pair(const char *, int); diff --git a/src/image.cpp b/src/image.cpp index f5be0bea4e..ab6824d7fb 100644 --- a/src/image.cpp +++ b/src/image.cpp @@ -1065,8 +1065,7 @@ void Image::write_PNG(FILE *fp) png_set_text(png_ptr,info_ptr,text_ptr,1); png_write_info(png_ptr,info_ptr); - png_bytep *row_pointers; - row_pointers = new png_bytep[height]; + png_bytep row_pointers[height]; for (int i=0; i < height; ++i) row_pointers[i] = (png_bytep) &writeBuffer[(height-i-1)*3*width]; @@ -1074,7 +1073,6 @@ void Image::write_PNG(FILE *fp) png_write_end(png_ptr, info_ptr); png_destroy_write_struct(&png_ptr, &info_ptr); - delete[] row_pointers; #endif } diff --git a/src/lattice.cpp b/src/lattice.cpp index d613b9e245..e5081d75f7 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -18,8 +18,8 @@ #include "update.h" #include "domain.h" #include "comm.h" -#include "memory.h" #include "force.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; diff --git a/src/memory.h b/src/memory.h index 9f4ccea23e..27b50f6198 100644 --- a/src/memory.h +++ b/src/memory.h @@ -25,7 +25,6 @@ namespace LAMMPS_NS { class Memory : protected Pointers { public: Memory(class LAMMPS *); - virtual ~Memory() {} void *smalloc(bigint n, const char *); void *srealloc(void *, bigint n, const char *); diff --git a/src/pair.cpp b/src/pair.cpp index 95a46ee441..a87b674608 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -223,10 +223,10 @@ void Pair::init() // I,I coeffs must be set // init_one() will check if I,J is set explicitly or inferred by mixing - if (!allocated) error->all(FLERR,"Not all pair coeffs are set"); + if (!allocated) error->all(FLERR,"All pair coeffs are not set"); for (i = 1; i <= atom->ntypes; i++) - if (setflag[i][i] == 0) error->all(FLERR,"Not all pair coeffs are set"); + if (setflag[i][i] == 0) error->all(FLERR,"All pair coeffs are not set"); // style-specific initialization diff --git a/src/pair.h b/src/pair.h index 35322e13f4..d1a41fbee6 100644 --- a/src/pair.h +++ b/src/pair.h @@ -277,10 +277,10 @@ This is likely not what you want to do. The exclusion settings will eliminate neighbors in the neighbor list, which the manybody potential needs to calculated its terms correctly. -E: Not all pair coeffs are set +E: All pair coeffs are not set -Pair coefficients must be set for all pairs of atom types in either the -data file or by the pair_coeff command before running a simulation. +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. E: Fix adapt interface to this pair style not supported diff --git a/src/random_mars.cpp b/src/random_mars.cpp index 8ce91127da..fb5725b8c3 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -17,7 +17,6 @@ #include "math.h" #include "random_mars.h" #include "error.h" -#include "math_inline.h" using namespace LAMMPS_NS; @@ -105,8 +104,7 @@ double RanMars::gaussian() rsq = v1*v1 + v2*v2; if (rsq < 1.0 && rsq != 0.0) again = 0; } - // fac = sqrt(-2.0*log(rsq)/rsq); - fac = MathInline::sqrtlgx2byx(rsq); + fac = sqrt(-2.0*log(rsq)/rsq); second = v1*fac; first = v2*fac; save = 1; diff --git a/src/random_park.cpp b/src/random_park.cpp index a5b3ed9106..89e89fc4b9 100644 --- a/src/random_park.cpp +++ b/src/random_park.cpp @@ -15,7 +15,6 @@ #include "math.h" #include "random_park.h" -#include "math_inline.h" #include "error.h" using namespace LAMMPS_NS; @@ -65,8 +64,7 @@ double RanPark::gaussian() rsq = v1*v1 + v2*v2; if (rsq < 1.0 && rsq != 0.0) again = 0; } - // fac = sqrt(-2.0*log(rsq)/rsq); - fac = MathInline::sqrtlgx2byx(rsq); + fac = sqrt(-2.0*log(rsq)/rsq); second = v1*fac; first = v2*fac; save = 1; diff --git a/src/universe.h b/src/universe.h index acda345ab9..897b6c9681 100644 --- a/src/universe.h +++ b/src/universe.h @@ -42,7 +42,7 @@ class Universe : protected Pointers { // proc uni2orig[I] in original communicator Universe(class LAMMPS *, MPI_Comm); - virtual ~Universe(); + ~Universe(); void reorder(char *, char *); void add_world(char *); int consistent(); diff --git a/src/update.h b/src/update.h index 2103395ef5..38bf1b7db9 100644 --- a/src/update.h +++ b/src/update.h @@ -47,7 +47,7 @@ class Update : protected Pointers { char *minimize_style; Update(class LAMMPS *); - virtual ~Update(); + ~Update(); void init(); void set_units(const char *); void create_integrate(int, char **, int);