diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 6ef6759f7f..d20b359fb2 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -671,7 +671,7 @@ void FixPour::pre_exchange() if (ninserted_atoms) { atom->natoms += ninserted_atoms; - if (atom->natoms < 0 || atom->natoms > MAXBIGINT) + if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); if (mode == MOLECULE) { atom->nbonds += onemols[imol]->nbonds * ninserted_mols; diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp index c107fc85c7..2f0874cf11 100644 --- a/src/KSPACE/ewald_disp.cpp +++ b/src/KSPACE/ewald_disp.cpp @@ -52,7 +52,7 @@ EwaldDisp::EwaldDisp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) ewaldflag = dispersionflag = dipoleflag = 1; accuracy_relative = fabs(force->numeric(FLERR,arg[0])); - memset(function, 0, EWALD_NORDER*sizeof(int)); + memset(function, 0, EWALD_NFUNCS*sizeof(int)); kenergy = kvirial = NULL; cek_local = cek_global = NULL; ekr_local = NULL; diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 7134ed4326..fcb483b680 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -245,7 +245,7 @@ void PPPMDisp::init() int *ptr = pair ? (int *) pair->extract("ewald_order",tmp) : NULL; double *p_cutoff = pair ? (double *) pair->extract("cut_coul",tmp) : NULL; double *p_cutoff_lj = pair ? (double *) pair->extract("cut_LJ",tmp) : NULL; - if (!(ptr||*p_cutoff||*p_cutoff_lj)) + if (!(ptr||p_cutoff||p_cutoff_lj)) error->all(FLERR,"KSpace style is incompatible with Pair style"); cutoff = *p_cutoff; cutoff_lj = *p_cutoff_lj; diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 2afd694dfd..58bbeffbfc 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -620,19 +620,19 @@ void PairAIREBO::FLJ(int eflag, int vflag) if (wik*wkj > best) { best = wik*wkj; npath = 3; - atomk = k; - delikS[0] = delik[0]; - delikS[1] = delik[1]; - delikS[2] = delik[2]; - rikS = rik; - wikS = wik; - dwikS = dwik; - deljkS[0] = deljk[0]; - deljkS[1] = deljk[1]; - deljkS[2] = deljk[2]; - rkjS = rkj; - wkjS = wkj; - dwkjS = dwkj; + atomk = k; + delikS[0] = delik[0]; + delikS[1] = delik[1]; + delikS[2] = delik[2]; + rikS = rik; + wikS = wik; + dwikS = dwik; + deljkS[0] = deljk[0]; + deljkS[1] = deljk[1]; + deljkS[2] = deljk[2]; + rkjS = rkj; + wkjS = wkj; + dwkjS = dwkj; if (best == 1.0) { done = 1; break; @@ -803,7 +803,7 @@ void PairAIREBO::FLJ(int eflag, int vflag) if (vflag_atom) v_tally3(atomi,atomj,atomk,fi,fj,delikS,deljkS); - } else { + } else if (npath == 4) { fpair1 = dC*dwikS*wkmS*wmjS / rikS; fi[0] = delikS[0]*fpair1; fi[1] = delikS[1]*fpair1; diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp index 94fbe62a99..20febf3a37 100644 --- a/src/MANYBODY/pair_bop.cpp +++ b/src/MANYBODY/pair_bop.cpp @@ -5108,11 +5108,7 @@ void PairBOP::read_table(char *filename) fgets(s,MAXLINE,fp); sscanf(s,"%lf%lf",&sigma_delta[i],&pi_delta[i]); fgets(s,MAXLINE,fp); - if(nws==3) { - sscanf(s,"%lf%lf%lf",&sigma_f[i],&sigma_k[i],&small3[i]); - } else { - sscanf(s,"%lf%lf%lf",&sigma_f[i],&sigma_k[i],&small3[i]); - } + sscanf(s,"%lf%lf%lf",&sigma_f[i],&sigma_k[i],&small3[i]); } if(nws==3) { for(i=0;idestroy(cutsq); delete [] map; delete [] type2frho; + map = NULL; + type2frho = NULL; memory->destroy(type2rhor); memory->destroy(type2z2r); } @@ -91,6 +93,7 @@ PairEAM::~PairEAM() memory->destroy(funcfl[i].zr); } memory->sfree(funcfl); + funcfl = NULL; } if (setfl) { @@ -101,6 +104,7 @@ PairEAM::~PairEAM() memory->destroy(setfl->rhor); memory->destroy(setfl->z2r); delete setfl; + setfl = NULL; } if (fs) { @@ -111,6 +115,7 @@ PairEAM::~PairEAM() memory->destroy(fs->rhor); memory->destroy(fs->z2r); delete fs; + fs = NULL; } memory->destroy(frho); diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index d65ea8bda0..d41e3ef325 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -1308,7 +1308,7 @@ void FixGCMC::attempt_molecule_insertion() fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); atom->natoms += natoms_per_molecule; - if (atom->natoms < 0 || atom->natoms > MAXBIGINT) + if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); atom->nbonds += onemols[imol]->nbonds; atom->nangles += onemols[imol]->nangles; @@ -1957,7 +1957,7 @@ void FixGCMC::attempt_molecule_insertion_full() fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); atom->natoms += natoms_per_molecule; - if (atom->natoms < 0 || atom->natoms > MAXBIGINT) + if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); atom->nbonds += onemols[imol]->nbonds; atom->nangles += onemols[imol]->nangles; diff --git a/src/MISC/fix_deposit.cpp b/src/MISC/fix_deposit.cpp index 0b35627ebe..7aa3f72ce0 100644 --- a/src/MISC/fix_deposit.cpp +++ b/src/MISC/fix_deposit.cpp @@ -563,7 +563,7 @@ void FixDeposit::pre_exchange() if (success) { atom->natoms += natom; - if (atom->natoms < 0 || atom->natoms > MAXBIGINT) + if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); if (mode == MOLECULE) { atom->nbonds += onemols[imol]->nbonds; diff --git a/src/MISC/fix_oneway.cpp b/src/MISC/fix_oneway.cpp index 2ed15f81ed..2f25b714b0 100644 --- a/src/MISC/fix_oneway.cpp +++ b/src/MISC/fix_oneway.cpp @@ -43,7 +43,7 @@ FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (nevery < 1) error->all(FLERR,"Illegal fix oneway command"); int len = strlen(arg[4]); - regionstr = new char[len]; + regionstr = new char[len+1]; strcpy(regionstr,arg[4]); if (strcmp(arg[5], "x") == 0) direction = X|PLUS; diff --git a/src/MOLECULE/angle_cosine_periodic.cpp b/src/MOLECULE/angle_cosine_periodic.cpp index b288c47029..8198f3d4b4 100644 --- a/src/MOLECULE/angle_cosine_periodic.cpp +++ b/src/MOLECULE/angle_cosine_periodic.cpp @@ -57,7 +57,7 @@ void AngleCosinePeriodic::compute(int eflag, int vflag) int i,i1,i2,i3,n,m,type,b_factor; double delx1,dely1,delz1,delx2,dely2,delz2; double eangle,f1[3],f3[3]; - double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22; + double rsq1,rsq2,r1,r2,c,a,a11,a12,a22; double tn,tn_1,tn_2,un,un_1,un_2; eangle = 0.0; @@ -120,10 +120,6 @@ void AngleCosinePeriodic::compute(int eflag, int vflag) un_1 = 2.0; un_2 = 0.0; - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - s = 1.0/s; - // force & energy tn_2 = c; diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp index 0eadf1686a..28b3fe8924 100644 --- a/src/MOLECULE/improper_umbrella.cpp +++ b/src/MOLECULE/improper_umbrella.cpp @@ -145,8 +145,8 @@ void ImproperUmbrella::compute(int eflag, int vflag) } } - if (c > 1.0) s = 1.0; - if (c < -1.0) s = -1.0; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/OPT/pair_eam_alloy_opt.h b/src/OPT/pair_eam_alloy_opt.h index c6056e3cff..3fdfc449e9 100644 --- a/src/OPT/pair_eam_alloy_opt.h +++ b/src/OPT/pair_eam_alloy_opt.h @@ -28,6 +28,7 @@ namespace LAMMPS_NS { class PairEAMAlloyOpt : public PairEAMAlloy, public PairEAMOpt { public: PairEAMAlloyOpt(class LAMMPS *); + virtual ~PairEAMAlloyOpt() {} }; } diff --git a/src/OPT/pair_eam_fs_opt.h b/src/OPT/pair_eam_fs_opt.h index f8fba98717..f0288d5c57 100644 --- a/src/OPT/pair_eam_fs_opt.h +++ b/src/OPT/pair_eam_fs_opt.h @@ -28,6 +28,7 @@ namespace LAMMPS_NS { class PairEAMFSOpt : public PairEAMFS, public PairEAMOpt { public: PairEAMFSOpt(class LAMMPS *); + virtual ~PairEAMFSOpt() {} }; } diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 8cb0be7462..6e5399ba82 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -207,7 +207,7 @@ void NEB::run() update->endstep = update->laststep = update->firststep + n1steps; update->nsteps = n1steps; update->max_eval = n1steps; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many timesteps for NEB"); update->minimize->setup(); @@ -275,7 +275,7 @@ void NEB::run() update->endstep = update->laststep = update->firststep + n2steps; update->nsteps = n2steps; update->max_eval = n2steps; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); update->minimize->init(); diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp index de3447db98..57aaf485f4 100644 --- a/src/REPLICA/prd.cpp +++ b/src/REPLICA/prd.cpp @@ -220,7 +220,7 @@ void PRD::command(int narg, char **arg) update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + nsteps; update->restrict_output = 1; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); lmp->init(); @@ -568,7 +568,7 @@ void PRD::quench() update->whichflag = 2; update->nsteps = maxiter; update->endstep = update->laststep = update->firststep + maxiter; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many iterations"); // full init works diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index 80381a7a3b..0f46867f99 100644 --- a/src/REPLICA/tad.cpp +++ b/src/REPLICA/tad.cpp @@ -199,7 +199,7 @@ void TAD::command(int narg, char **arg) update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + nsteps; update->restrict_output = 1; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); lmp->init(); @@ -406,7 +406,7 @@ void TAD::command(int narg, char **arg) timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms); } - if (me_universe == 0) fclose(ulogfile_neb); + if ((me_universe == 0) && ulogfile_neb) fclose(ulogfile_neb); if (me == 0) { if (screen) fprintf(screen,"\nTAD done\n"); @@ -478,7 +478,7 @@ void TAD::quench() update->whichflag = 2; update->nsteps = maxiter; update->endstep = update->laststep = update->firststep + maxiter; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many iterations"); // full init works diff --git a/src/REPLICA/temper.cpp b/src/REPLICA/temper.cpp index 45e100ad52..8771ab5375 100644 --- a/src/REPLICA/temper.cpp +++ b/src/REPLICA/temper.cpp @@ -107,7 +107,7 @@ void Temper::command(int narg, char **arg) update->nsteps = nsteps; update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + nsteps; - if (update->laststep < 0 || update->laststep > MAXBIGINT) + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); lmp->init(); diff --git a/src/SHOCK/fix_append_atoms.cpp b/src/SHOCK/fix_append_atoms.cpp index 571e77b129..f01fbdfbbb 100644 --- a/src/SHOCK/fix_append_atoms.cpp +++ b/src/SHOCK/fix_append_atoms.cpp @@ -506,7 +506,7 @@ void FixAppendAtoms::pre_exchange() if (addtotal) { domain->reset_box(); atom->natoms += addtotal; - if (atom->natoms < 0 || atom->natoms > MAXBIGINT) + if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); if (atom->tag_enable) atom->tag_extend(); if (atom->map_style) { diff --git a/src/SHOCK/fix_wall_piston.cpp b/src/SHOCK/fix_wall_piston.cpp index 2f908bccb4..9217a24f6f 100644 --- a/src/SHOCK/fix_wall_piston.cpp +++ b/src/SHOCK/fix_wall_piston.cpp @@ -41,6 +41,7 @@ FixWallPiston::FixWallPiston(LAMMPS *lmp, int narg, char **arg) : if (narg < 4) error->all(FLERR,"Illegal fix wall/piston command"); randomt = NULL; + gfactor1 = gfactor2 = NULL; tempflag = 0; scaleflag = 1; roughflag = 0; @@ -150,6 +151,15 @@ FixWallPiston::FixWallPiston(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ +FixWallPiston::~FixWallPiston() +{ + delete[] gfactor2; + delete[] gfactor1; + delete randomt; +} + +/* ---------------------------------------------------------------------- */ + int FixWallPiston::setmask() { int mask = 0; diff --git a/src/SHOCK/fix_wall_piston.h b/src/SHOCK/fix_wall_piston.h index 8ad3c25b57..28135b2359 100644 --- a/src/SHOCK/fix_wall_piston.h +++ b/src/SHOCK/fix_wall_piston.h @@ -26,6 +26,7 @@ namespace LAMMPS_NS { class FixWallPiston : public Fix { public: FixWallPiston(class LAMMPS *, int, char **); + virtual ~FixWallPiston(); int setmask(); void post_integrate(); void initial_integrate(int); diff --git a/src/USER-MISC/compute_basal_atom.cpp b/src/USER-MISC/compute_basal_atom.cpp index 85a08d1838..cf90b50d8a 100644 --- a/src/USER-MISC/compute_basal_atom.cpp +++ b/src/USER-MISC/compute_basal_atom.cpp @@ -243,7 +243,7 @@ void ComputeBasalAtom::compute_peratom() j1[1]=2; j1[2]=2; } - xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0; + xmean5 = ymean5 = zmean5 = xmean6 = ymean6 = zmean6 = xmean7 = ymean7 = zmean7 = 0.0; for (j = 0; j < chi[0]; j++) { for (k = j+1; k < chi[0]; k++) { //get cross products @@ -261,27 +261,29 @@ void ComputeBasalAtom::compute_peratom() y7[count] = y4[count]*copysign(1.0,z4[count]); z7[count] = z4[count]*copysign(1.0,z4[count]); //get average cross products - xmean5 = xmean5 + x5[count]; - ymean5 = ymean5 + y5[count]; - zmean5 = zmean5 + z5[count]; - xmean6 = xmean6 + x6[count]; - ymean6 = ymean6 + y6[count]; - zmean6 = zmean6 + z6[count]; - xmean7 = xmean7 + x7[count]; - ymean7 = ymean7 + y7[count]; - zmean6 = zmean6 + z7[count]; + xmean5 += x5[count]; + ymean5 += y5[count]; + zmean5 += z5[count]; + xmean6 += x6[count]; + ymean6 += y6[count]; + zmean6 += z6[count]; + xmean7 += x7[count]; + ymean7 += y7[count]; + zmean6 += z7[count]; count++; } } - xmean5 = xmean5/count; - xmean6 = xmean6/count; - xmean7 = xmean7/count; - ymean5 = ymean5/count; - ymean6 = ymean6/count; - ymean7 = ymean7/count; - zmean5 = zmean5/count; - zmean6 = zmean6/count; - zmean7 = zmean7/count; + if (count > 0) { + xmean5 /= count; + xmean6 /= count; + xmean7 /= count; + ymean5 /= count; + ymean6 /= count; + ymean7 /= count; + zmean5 /= count; + zmean6 /= count; + zmean7 /= count; + } var5 = var6 = var7 = 0.0; //find standard deviations for (j=0;jffile[I] = ffile_tmp[i]; I++; } + + // clean up temporary storage + delete[] phifile_tmp; + delete[] ffile_tmp; + delete[] efile_tmp; } // spline read-in and compute r,e,f vectors within table @@ -1147,10 +1164,14 @@ void DihedralTable::spline_table(Table *tb) memory->create(tb->e2file,tb->ninput,"dihedral:e2file"); memory->create(tb->f2file,tb->ninput,"dihedral:f2file"); - cyc_spline(tb->phifile, tb->efile, tb->ninput, MY_2PI, tb->e2file); + if (cyc_spline(tb->phifile, tb->efile, tb->ninput, + MY_2PI,tb->e2file,comm->me == 0)) + error->one(FLERR,"Error computing dihedral spline tables"); if (! tb->f_unspecified) { - cyc_spline(tb->phifile, tb->ffile, tb->ninput, MY_2PI, tb->f2file); + if (cyc_spline(tb->phifile, tb->ffile, tb->ninput, + MY_2PI, tb->f2file, comm->me == 0)) + error->one(FLERR,"Error computing dihedral spline tables"); } // CHECK to help make sure the user calculated forces in a way @@ -1302,9 +1323,9 @@ void DihedralTable::compute_table(Table *tb) - cyc_spline(tb->phi, tb->e, tablength, MY_2PI, tb->e2); + cyc_spline(tb->phi, tb->e, tablength, MY_2PI, tb->e2, comm->me == 0); if (! tb->f_unspecified) - cyc_spline(tb->phi, tb->f, tablength, MY_2PI, tb->f2); + cyc_spline(tb->phi, tb->f, tablength, MY_2PI, tb->f2, comm->me == 0); } @@ -1339,13 +1360,13 @@ void DihedralTable::param_extract(Table *tb, char *line) else if (strcmp(word,"CHECKU") == 0) { word = strtok(NULL," \t\n\r\f"); memory->sfree(checkU_fname); - memory->create(checkU_fname,strlen(word),"dihedral_table:checkU"); + memory->create(checkU_fname,strlen(word)+1,"dihedral_table:checkU"); strcpy(checkU_fname, word); } else if (strcmp(word,"CHECKF") == 0) { word = strtok(NULL," \t\n\r\f"); memory->sfree(checkF_fname); - memory->create(checkF_fname,strlen(word),"dihedral_table:checkF"); + memory->create(checkF_fname,strlen(word)+1,"dihedral_table:checkF"); strcpy(checkF_fname, word); } // COMMENTING OUT: equilibrium angles are not supported diff --git a/src/USER-MISC/fix_ave_spatial_sphere.cpp b/src/USER-MISC/fix_ave_spatial_sphere.cpp index c654db7ef8..bf82173be2 100644 --- a/src/USER-MISC/fix_ave_spatial_sphere.cpp +++ b/src/USER-MISC/fix_ave_spatial_sphere.cpp @@ -827,7 +827,7 @@ void FixAveSpatialSphere::end_of_step() fflush(fp); if (overwrite) { long fileend = ftell(fp); - ftruncate(fileno(fp),fileend); + if (fileend > 0) ftruncate(fileno(fp),fileend); } } } diff --git a/src/USER-MISC/fix_ipi.cpp b/src/USER-MISC/fix_ipi.cpp index 58251a814e..07dbeb2500 100644 --- a/src/USER-MISC/fix_ipi.cpp +++ b/src/USER-MISC/fix_ipi.cpp @@ -18,6 +18,7 @@ #include "mpi.h" #include "stdio.h" #include "string.h" +#include "stdlib.h" #include "fix_ipi.h" #include "atom.h" #include "force.h" @@ -187,8 +188,8 @@ FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[1],"all")) error->warning(FLERR,"Fix ipi always uses group all"); - strcpy(host, arg[3]); - port=force->inumeric(FLERR,arg[4]); + host = strdup(arg[3]); + port = force->inumeric(FLERR,arg[4]); inet = ((narg > 5) && (strcmp(arg[5],"unix") ==0) ) ? 0 : 1; master = (comm->me==0) ? 1 : 0; @@ -218,6 +219,7 @@ FixIPI::FixIPI(LAMMPS *lmp, int narg, char **arg) : FixIPI::~FixIPI() { if (bsize) delete[] buffer; + free(host); modify->delete_compute("IPI_TEMP"); modify->delete_compute("IPI_PRESS"); } diff --git a/src/USER-MISC/fix_ipi.h b/src/USER-MISC/fix_ipi.h index 010f0ba17d..0859a513c5 100644 --- a/src/USER-MISC/fix_ipi.h +++ b/src/USER-MISC/fix_ipi.h @@ -34,7 +34,7 @@ class FixIPI : public Fix { virtual void final_integrate(); protected: - char host[1024]; int port; int inet, master, hasdata; + char *host; int port; int inet, master, hasdata; int ipisock, me; double *buffer; long bsize; int kspace_flag; }; diff --git a/src/USER-MISC/fix_srp.cpp b/src/USER-MISC/fix_srp.cpp index 833dc41d28..613204a55e 100644 --- a/src/USER-MISC/fix_srp.cpp +++ b/src/USER-MISC/fix_srp.cpp @@ -178,7 +178,9 @@ void FixSRP::setup_pre_force(int zz) for (int n = 0; n < nbondlist; n++) { // consider only the user defined bond type - if(bondlist[n][2] != btype) continue; + // btype of zero considers all bonds + if(btype > 0 && bondlist[n][2] != btype) + continue; i = bondlist[n][0]; j = bondlist[n][1]; diff --git a/src/USER-MISC/fix_ti_spring.cpp b/src/USER-MISC/fix_ti_spring.cpp index 139795557c..529e650978 100755 --- a/src/USER-MISC/fix_ti_spring.cpp +++ b/src/USER-MISC/fix_ti_spring.cpp @@ -76,8 +76,8 @@ FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) : t_switch = atoi(arg[4]); // Switching time. t_equil = atoi(arg[5]); // Equilibration time. t0 = update->ntimestep; // Initial time. - if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/spring command"); - if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/spring command"); + if (t_switch <= 0.0) error->all(FLERR,"Illegal fix ti/spring command"); + if (t_equil <= 0.0) error->all(FLERR,"Illegal fix ti/spring command"); // Coupling parameter initialization. sf = 1; diff --git a/src/USER-MISC/fix_ttm_mod.cpp b/src/USER-MISC/fix_ttm_mod.cpp index d66e49d6a2..3a6aafe700 100644 --- a/src/USER-MISC/fix_ttm_mod.cpp +++ b/src/USER-MISC/fix_ttm_mod.cpp @@ -225,6 +225,10 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) : t_surface_l = surface_l; mult_factor = intensity; duration = 0.0; + v_0_sq = v_0*v_0; + // error checks + if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0) + error->all(FLERR,"Fix ttm number of nodes must be > 0"); surface_double = double(t_surface_l)*(domain->xprd/nxnodes); if ((C_limit+esheat_0) < 0.0) error->all(FLERR,"Fix ttm electronic_specific_heat must be >= 0.0"); @@ -234,10 +238,6 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) : if (gamma_s < 0.0) error->all(FLERR,"Fix ttm gamma_s must be >= 0.0"); if (v_0 < 0.0) error->all(FLERR,"Fix ttm v_0 must be >= 0.0"); if (ionic_density <= 0.0) error->all(FLERR,"Fix ttm ionic_density must be > 0.0"); - v_0_sq = v_0*v_0; - // error check - if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0) - error->all(FLERR,"Fix ttm number of nodes must be > 0"); if (seed <= 0) error->all(FLERR,"Invalid random number seed in fix ttm command"); if (surface_l < 0) error->all(FLERR,"Surface coordinates must be >= 0"); if (surface_l >= surface_r) error->all(FLERR, "Left surface coordinate must be less than right surface coordinate"); diff --git a/src/USER-MISC/improper_fourier.cpp b/src/USER-MISC/improper_fourier.cpp index edfd95f20f..f3fe065d6e 100644 --- a/src/USER-MISC/improper_fourier.cpp +++ b/src/USER-MISC/improper_fourier.cpp @@ -173,8 +173,8 @@ void ImproperFourier::addone(const int &i1,const int &i2,const int &i3,const int } } - if (c > 1.0) s = 1.0; - if (c < -1.0) s = -1.0; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/USER-MISC/pair_srp.cpp b/src/USER-MISC/pair_srp.cpp index ec7730a1f2..fedbe0a302 100644 --- a/src/USER-MISC/pair_srp.cpp +++ b/src/USER-MISC/pair_srp.cpp @@ -59,6 +59,8 @@ static const char cite_srp[] = " pages = {134903}\n" "}\n\n"; +static int srp_instance = 0; + /* ---------------------------------------------------------------------- set size of pair comms in constructor ---------------------------------------------------------------------- */ @@ -71,6 +73,10 @@ PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) nextra = 1; segment = NULL; + fix_id = strdup("XX_FIX_SRP"); + fix_id[0] = '0' + srp_instance / 10; + fix_id[1] = '0' + srp_instance % 10; + ++srp_instance; } /* ---------------------------------------------------------------------- @@ -112,9 +118,10 @@ PairSRP::~PairSRP() } // check nfix in case all fixes have already been deleted - if (modify->nfix) modify->delete_fix("mysrp"); -} - + if (modify->nfix) modify->delete_fix(fix_id); + free(fix_id); +} + /* ---------------------------------------------------------------------- compute bond-bond repulsions ------------------------------------------------------------------------- */ @@ -336,8 +343,14 @@ void PairSRP::settings(int narg, char **arg) error->all(FLERR,"Illegal pair_style command"); cut_global = force->numeric(FLERR,arg[0]); - btype = force->inumeric(FLERR,arg[1]); - if (btype > atom->nbondtypes) error->all(FLERR,"Illegal pair_style command"); + // wildcard + if (strcmp(arg[1],"*") == 0) + btype = 0; + else { + btype = force->inumeric(FLERR,arg[1]); + if ((btype > atom->nbondtypes) || (btype <= 0)) + error->all(FLERR,"Illegal pair_style command"); + } // settings midpoint = 0; @@ -426,7 +439,7 @@ void PairSRP::init_style() // invoke here instead of constructor // to make restart possible char **fixarg = new char*[3]; - fixarg[0] = (char *) "mysrp"; + fixarg[0] = fix_id; fixarg[1] = (char *) "all"; fixarg[2] = (char *) "SRP"; modify->add_fix(3,fixarg); diff --git a/src/USER-MISC/pair_srp.h b/src/USER-MISC/pair_srp.h index 998f8a17e2..ecd352cdeb 100644 --- a/src/USER-MISC/pair_srp.h +++ b/src/USER-MISC/pair_srp.h @@ -54,6 +54,7 @@ class PairSRP : public Pair { int bptype; int btype; class Fix *f_srp; + char *fix_id; int exclude,maxcount; int **segment; }; diff --git a/src/USER-MOLFILE/molfile_interface.cpp b/src/USER-MOLFILE/molfile_interface.cpp index d0aff40f32..b9ee2238a9 100644 --- a/src/USER-MOLFILE/molfile_interface.cpp +++ b/src/USER-MOLFILE/molfile_interface.cpp @@ -756,10 +756,13 @@ int MolfileInterface::timestep(float *coords, float *vels, *simtime = t->physical_time; } - if (rv == MOLFILE_EOF) + if (rv == MOLFILE_EOF) { + delete t; return 1; + } } + delete t; return 0; } diff --git a/src/USER-OMP/angle_cosine_periodic_omp.cpp b/src/USER-OMP/angle_cosine_periodic_omp.cpp index 311152f06c..3bd4ef1212 100644 --- a/src/USER-OMP/angle_cosine_periodic_omp.cpp +++ b/src/USER-OMP/angle_cosine_periodic_omp.cpp @@ -91,7 +91,7 @@ void AngleCosinePeriodicOMP::eval(int nfrom, int nto, ThrData * const thr) int i,i1,i2,i3,n,m,type,b_factor; double delx1,dely1,delz1,delx2,dely2,delz2; double eangle,f1[3],f3[3]; - double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22; + double rsq1,rsq2,r1,r2,c,a,a11,a12,a22; double tn,tn_1,tn_2,un,un_1,un_2; const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; @@ -149,10 +149,6 @@ void AngleCosinePeriodicOMP::eval(int nfrom, int nto, ThrData * const thr) un_1 = 2.0; un_2 = 0.0; - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - s = 1.0/s; - // force & energy tn_2 = c; diff --git a/src/USER-OMP/improper_fourier_omp.cpp b/src/USER-OMP/improper_fourier_omp.cpp index 49fcef23d4..6c32c8e538 100644 --- a/src/USER-OMP/improper_fourier_omp.cpp +++ b/src/USER-OMP/improper_fourier_omp.cpp @@ -206,8 +206,8 @@ void ImproperFourierOMP::add1_thr(const int i1,const int i2, } } - if (c > 1.0) s = 1.0; - if (c < -1.0) s = -1.0; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/USER-OMP/improper_umbrella_omp.cpp b/src/USER-OMP/improper_umbrella_omp.cpp index 689dbdfe65..f100d3a0a9 100644 --- a/src/USER-OMP/improper_umbrella_omp.cpp +++ b/src/USER-OMP/improper_umbrella_omp.cpp @@ -171,8 +171,8 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr) } } - if (c > 1.0) s = 1.0; - if (c < -1.0) s = -1.0; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp index ee36e5610e..e61611d41c 100644 --- a/src/USER-OMP/neigh_full_omp.cpp +++ b/src/USER-OMP/neigh_full_omp.cpp @@ -599,7 +599,6 @@ void Neighbor::full_multi_omp(NeighList *list) onemols[imol]->nspecial[iatom], tag[j]-tagprev); else which = 0; - which = find_special(special[i],nspecial[i],tag[j]); if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index c63ac3cd35..f2149088c1 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -2067,19 +2067,19 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, if (wik*wkj > best) { best = wik*wkj; npath = 3; - atomk = k; - delikS[0] = delik[0]; - delikS[1] = delik[1]; - delikS[2] = delik[2]; - rikS = rik; - wikS = wik; - dwikS = dwik; - deljkS[0] = deljk[0]; - deljkS[1] = deljk[1]; - deljkS[2] = deljk[2]; - rkjS = rkj; - wkjS = wkj; - dwkjS = dwkj; + atomk = k; + delikS[0] = delik[0]; + delikS[1] = delik[1]; + delikS[2] = delik[2]; + rikS = rik; + wikS = wik; + dwikS = dwik; + deljkS[0] = deljk[0]; + deljkS[1] = deljk[1]; + deljkS[2] = deljk[2]; + rkjS = rkj; + wkjS = wkj; + dwkjS = dwkj; if (best == 1.0) { done = 1; break; @@ -2128,13 +2128,13 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, atomm = m; delkmS[0] = delkm[0]; delkmS[1] = delkm[1]; - delkmS[2] = delkm[2]; + delkmS[2] = delkm[2]; rkmS = rkm; wkmS = wkm; dwkmS = dwkm; deljmS[0] = deljm[0]; deljmS[1] = deljm[1]; - deljmS[2] = deljm[2]; + deljmS[2] = deljm[2]; rmjS = rmj; wmjS = wmj; dwmjS = dwmj; @@ -2250,7 +2250,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, if (vflag_atom) v_tally3_thr(atomi,atomj,atomk,fi,fj,delikS,deljkS,thr); - } else { + } else if (npath == 4) { fpair1 = dC*dwikS*wkmS*wmjS / rikS; fi[0] = delikS[0]*fpair1; fi[1] = delikS[1]*fpair1; diff --git a/src/USER-OMP/pair_brownian_omp.cpp b/src/USER-OMP/pair_brownian_omp.cpp index 56d9dee308..8412644ff3 100644 --- a/src/USER-OMP/pair_brownian_omp.cpp +++ b/src/USER-OMP/pair_brownian_omp.cpp @@ -48,6 +48,7 @@ PairBrownianOMP::PairBrownianOMP(LAMMPS *lmp) : suffix_flag |= Suffix::OMP; respa_enable = 0; random_thr = NULL; + nthreads = 0; } /* ---------------------------------------------------------------------- */ @@ -55,7 +56,7 @@ PairBrownianOMP::PairBrownianOMP(LAMMPS *lmp) : PairBrownianOMP::~PairBrownianOMP() { if (random_thr) { - for (int i=1; i < comm->nthreads; ++i) + for (int i=1; i < nthreads; ++i) delete random_thr[i]; delete[] random_thr; @@ -72,7 +73,6 @@ void PairBrownianOMP::compute(int eflag, int vflag) } else evflag = vflag_fdotr = 0; const int nall = atom->nlocal + atom->nghost; - const int nthreads = comm->nthreads; const int inum = list->inum; // This section of code adjusts R0/RT0/RS0 if necessary due to changes @@ -116,13 +116,24 @@ void PairBrownianOMP::compute(int eflag, int vflag) } } + // number of threads has changed. reallocate pool of pRNGs + if (nthreads != comm->nthreads) { + if (random_thr) { + for (int i=1; i < nthreads; ++i) + delete random_thr[i]; - if (!random_thr) + delete[] random_thr; + } + + nthreads = comm->nthreads; random_thr = new RanMars*[nthreads]; + for (int i=1; i < nthreads; ++i) + random_thr[i] = NULL; - // to ensure full compatibility with the serial Brownian style - // we use is random number generator instance for thread 0 - random_thr[0] = random; + // to ensure full compatibility with the serial Brownian style + // we use is random number generator instance for thread 0 + random_thr[0] = random; + } #if defined(_OPENMP) #pragma omp parallel default(none) shared(eflag,vflag) @@ -137,9 +148,10 @@ void PairBrownianOMP::compute(int eflag, int vflag) // generate a random number generator instance for // all threads != 0. make sure we use unique seeds. - if (random_thr && tid > 0) + if ((tid > 0) && (random_thr[tid] == NULL)) random_thr[tid] = new RanMars(Pair::lmp, seed + comm->me + comm->nprocs*tid); + if (flaglog) { if (evflag) { if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); @@ -395,8 +407,8 @@ double PairBrownianOMP::memory_usage() { double bytes = memory_usage_thr(); bytes += PairBrownian::memory_usage(); - bytes += comm->nthreads * sizeof(RanMars*); - bytes += comm->nthreads * sizeof(RanMars); + bytes += nthreads * sizeof(RanMars*); + bytes += nthreads * sizeof(RanMars); return bytes; } diff --git a/src/USER-OMP/pair_brownian_omp.h b/src/USER-OMP/pair_brownian_omp.h index cdabe08ba1..c49560cbca 100644 --- a/src/USER-OMP/pair_brownian_omp.h +++ b/src/USER-OMP/pair_brownian_omp.h @@ -40,6 +40,7 @@ class PairBrownianOMP : public PairBrownian, public ThrOMP { protected: class RanMars **random_thr; + int nthreads; private: template diff --git a/src/USER-OMP/pair_brownian_poly_omp.cpp b/src/USER-OMP/pair_brownian_poly_omp.cpp index 0ac97a935b..e0124204de 100644 --- a/src/USER-OMP/pair_brownian_poly_omp.cpp +++ b/src/USER-OMP/pair_brownian_poly_omp.cpp @@ -48,6 +48,7 @@ PairBrownianPolyOMP::PairBrownianPolyOMP(LAMMPS *lmp) : suffix_flag |= Suffix::OMP; respa_enable = 0; random_thr = NULL; + nthreads = 0; } /* ---------------------------------------------------------------------- */ @@ -55,7 +56,7 @@ PairBrownianPolyOMP::PairBrownianPolyOMP(LAMMPS *lmp) : PairBrownianPolyOMP::~PairBrownianPolyOMP() { if (random_thr) { - for (int i=1; i < comm->nthreads; ++i) + for (int i=1; i < nthreads; ++i) delete random_thr[i]; delete[] random_thr; @@ -72,7 +73,6 @@ void PairBrownianPolyOMP::compute(int eflag, int vflag) } else evflag = vflag_fdotr = 0; const int nall = atom->nlocal + atom->nghost; - const int nthreads = comm->nthreads; const int inum = list->inum; // This section of code adjusts R0/RT0/RS0 if necessary due to changes @@ -117,12 +117,24 @@ void PairBrownianPolyOMP::compute(int eflag, int vflag) } - if (!random_thr) - random_thr = new RanMars*[nthreads]; + // number of threads has changed. reallocate pool of pRNGs + if (nthreads != comm->nthreads) { + if (random_thr) { + for (int i=1; i < nthreads; ++i) + delete random_thr[i]; - // to ensure full compatibility with the serial BrownianPoly style - // we use is random number generator instance for thread 0 - random_thr[0] = random; + delete[] random_thr; + } + + nthreads = comm->nthreads; + random_thr = new RanMars*[nthreads]; + for (int i=1; i < nthreads; ++i) + random_thr[i] = NULL; + + // to ensure full compatibility with the serial BrownianPoly style + // we use is random number generator instance for thread 0 + random_thr[0] = random; + } #if defined(_OPENMP) #pragma omp parallel default(none) shared(eflag,vflag) @@ -137,9 +149,10 @@ void PairBrownianPolyOMP::compute(int eflag, int vflag) // generate a random number generator instance for // all threads != 0. make sure we use unique seeds. - if (random_thr && tid > 0) + if ((tid > 0) && (random_thr[tid] == NULL)) random_thr[tid] = new RanMars(Pair::lmp, seed + comm->me + comm->nprocs*tid); + if (flaglog) { if (evflag) eval<1,1>(ifrom, ito, thr); @@ -384,8 +397,8 @@ double PairBrownianPolyOMP::memory_usage() { double bytes = memory_usage_thr(); bytes += PairBrownianPoly::memory_usage(); - bytes += comm->nthreads * sizeof(RanMars*); - bytes += comm->nthreads * sizeof(RanMars); + bytes += nthreads * sizeof(RanMars*); + bytes += nthreads * sizeof(RanMars); return bytes; } diff --git a/src/USER-OMP/pair_brownian_poly_omp.h b/src/USER-OMP/pair_brownian_poly_omp.h index 663696c6c6..df2c0b34bd 100644 --- a/src/USER-OMP/pair_brownian_poly_omp.h +++ b/src/USER-OMP/pair_brownian_poly_omp.h @@ -40,6 +40,7 @@ class PairBrownianPolyOMP : public PairBrownianPoly, public ThrOMP { protected: class RanMars **random_thr; + int nthreads; private: template diff --git a/src/USER-OMP/pair_dpd_omp.cpp b/src/USER-OMP/pair_dpd_omp.cpp index 23408469e3..625d5bc19f 100644 --- a/src/USER-OMP/pair_dpd_omp.cpp +++ b/src/USER-OMP/pair_dpd_omp.cpp @@ -74,12 +74,13 @@ void PairDPDOMP::compute(int eflag, int vflag) nthreads = comm->nthreads; random_thr = new RanMars*[nthreads]; for (int i=1; i < nthreads; ++i) - random_thr[i] = NULL; + random_thr[i] = NULL; // to ensure full compatibility with the serial DPD style // we use the serial random number generator instance for thread 0 random_thr[0] = random; } + #if defined(_OPENMP) #pragma omp parallel default(none) shared(eflag,vflag) #endif @@ -93,7 +94,7 @@ void PairDPDOMP::compute(int eflag, int vflag) // generate a random number generator instance for // all threads != 0. make sure we use unique seeds. - if (random_thr && (tid > 0) && (random_thr[tid] == NULL)) + if ((tid > 0) && (random_thr[tid] == NULL)) random_thr[tid] = new RanMars(Pair::lmp, seed + comm->me + comm->nprocs*tid); diff --git a/src/USER-OMP/pair_dpd_tstat_omp.cpp b/src/USER-OMP/pair_dpd_tstat_omp.cpp index ca5b70186b..37b2994e3d 100644 --- a/src/USER-OMP/pair_dpd_tstat_omp.cpp +++ b/src/USER-OMP/pair_dpd_tstat_omp.cpp @@ -93,7 +93,7 @@ void PairDPDTstatOMP::compute(int eflag, int vflag) // generate a random number generator instance for // all threads != 0. make sure we use unique seeds. - if (random_thr && (tid > 0) && (random_thr[tid] == NULL)) + if ((tid > 0) && (random_thr[tid] == NULL)) random_thr[tid] = new RanMars(Pair::lmp, seed + comm->me + comm->nprocs*tid); diff --git a/src/USER-OMP/thr_data.h b/src/USER-OMP/thr_data.h index 3f1d866a80..3cea1b0018 100644 --- a/src/USER-OMP/thr_data.h +++ b/src/USER-OMP/thr_data.h @@ -35,7 +35,7 @@ class ThrData { public: ThrData(int tid, class Timer *t); - ~ThrData() {}; + ~ThrData() { delete _timer; _timer = NULL; }; void check_tid(int); // thread id consistency check int get_tid() const { return _tid; }; // our thread id. @@ -140,7 +140,7 @@ class ThrData { // disabled default methods private: - ThrData() : _tid(-1) {}; + ThrData() : _tid(-1), _timer(NULL) {}; }; //////////////////////////////////////////////////////////////////////// diff --git a/src/USER-PHONON/fix_phonon.h b/src/USER-PHONON/fix_phonon.h index 891235cfec..f1a8af4812 100644 --- a/src/USER-PHONON/fix_phonon.h +++ b/src/USER-PHONON/fix_phonon.h @@ -31,6 +31,17 @@ FixStyle(phonon,FixPhonon) #ifndef FIX_PHONON_H #define FIX_PHONON_H +#include "lmptype.h" +#include "mpi.h" + +#ifdef FFT_SINGLE +typedef float FFT_SCALAR; +#define MPI_FFT_SCALAR MPI_FLOAT +#else +typedef double FFT_SCALAR; +#define MPI_FFT_SCALAR MPI_DOUBLE +#endif + #include #include "fix.h" #include @@ -72,7 +83,7 @@ class FixPhonon : public Fix { int mynpt,mynq,fft_nsend; int *fft_cnts, *fft_disp; int fft_dim, fft_dim2; - double *fft_data; + FFT_SCALAR *fft_data; tagint itag; // index variables int idx, idq; // more index variables diff --git a/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp b/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp index 5c20b803e8..79c22e1056 100644 --- a/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp +++ b/src/USER-SMD/fix_smd_tlsph_reference_configuration.cpp @@ -330,6 +330,7 @@ void FixSMD_TLSPH_ReferenceConfiguration::setup(int vflag) { int nall, countall; MPI_Allreduce(&n, &nall, 1, MPI_INT, MPI_SUM, world); MPI_Allreduce(&count, &countall, 1, MPI_INT, MPI_SUM, world); + if (countall < 1) countall = 1; if (comm->me == 0) { if (screen) { diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index 7c5f433bb5..65bd4eabe2 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -95,6 +95,7 @@ ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) : typecount = new int[ntypes+1]; icount = new int[npairs]; jcount = new int[npairs]; + duplicates = new int[npairs]; } /* ---------------------------------------------------------------------- */ @@ -113,13 +114,14 @@ ComputeRDF::~ComputeRDF() delete [] typecount; delete [] icount; delete [] jcount; + delete [] duplicates; } /* ---------------------------------------------------------------------- */ void ComputeRDF::init() { - int i,m; + int i,j,m; if (force->pair) delr = force->pair->cutforce / nbin; else error->all(FLERR,"Compute rdf requires a pair style be defined"); @@ -143,12 +145,17 @@ void ComputeRDF::init() // icount = # of I atoms participating in I,J pairs for each histogram // jcount = # of J atoms participating in I,J pairs for each histogram + // duplicates = # of atoms in both groups I and J for each histogram for (m = 0; m < npairs; m++) { icount[m] = 0; for (i = ilo[m]; i <= ihi[m]; i++) icount[m] += typecount[i]; jcount[m] = 0; for (i = jlo[m]; i <= jhi[m]; i++) jcount[m] += typecount[i]; + duplicates[m] = 0; + for (i = ilo[m]; i <= ihi[m]; i++) + for (j = jlo[m]; j <= jhi[m]; j++) + if (i == j) duplicates[m] += typecount[i]; } int *scratch = new int[npairs]; @@ -156,6 +163,8 @@ void ComputeRDF::init() for (i = 0; i < npairs; i++) icount[i] = scratch[i]; MPI_Allreduce(jcount,scratch,npairs,MPI_INT,MPI_SUM,world); for (i = 0; i < npairs; i++) jcount[i] = scratch[i]; + MPI_Allreduce(duplicates,scratch,npairs,MPI_INT,MPI_SUM,world); + for (i = 0; i < npairs; i++) duplicates[i] = scratch[i]; delete [] scratch; // need an occasional half neighbor list @@ -265,25 +274,28 @@ void ComputeRDF::compute_array() MPI_Allreduce(hist[0],histall[0],npairs*nbin,MPI_DOUBLE,MPI_SUM,world); // convert counts to g(r) and coord(r) and copy into output array - // nideal = # of J atoms surrounding single I atom in a single bin - // assuming J atoms are at uniform density + // vfrac = fraction of volume in shell m + // npairs = number of pairs, corrected for duplicates + // duplicates = pairs in which both atoms are the same - double constant,nideal,gr,ncoord,rlower,rupper; + double constant,vfrac,gr,ncoord,rlower,rupper,normfac; if (domain->dimension == 3) { constant = 4.0*MY_PI / (3.0*domain->xprd*domain->yprd*domain->zprd); for (m = 0; m < npairs; m++) { + normfac = (icount[m] > 0) ? static_cast(jcount[m]) + - static_cast(duplicates[m])/icount[m] : 0.0; ncoord = 0.0; for (ibin = 0; ibin < nbin; ibin++) { rlower = ibin*delr; rupper = (ibin+1)*delr; - nideal = constant * - (rupper*rupper*rupper - rlower*rlower*rlower) * jcount[m]; - if (icount[m]*nideal != 0.0) - gr = histall[m][ibin] / (icount[m]*nideal); + vfrac = constant * (rupper*rupper*rupper - rlower*rlower*rlower); + if (vfrac * normfac != 0.0) + gr = histall[m][ibin] / (vfrac * normfac * icount[m]); else gr = 0.0; - ncoord += gr*nideal; + if (icount[m] != 0) + ncoord += gr * vfrac * normfac; array[ibin][1+2*m] = gr; array[ibin][2+2*m] = ncoord; } @@ -294,14 +306,17 @@ void ComputeRDF::compute_array() for (m = 0; m < npairs; m++) { ncoord = 0.0; + normfac = (icount[m] > 0) ? static_cast(jcount[m]) + - static_cast(duplicates[m])/icount[m] : 0.0; for (ibin = 0; ibin < nbin; ibin++) { rlower = ibin*delr; rupper = (ibin+1)*delr; - nideal = constant * (rupper*rupper - rlower*rlower) * jcount[m]; - if (icount[m]*nideal != 0.0) - gr = histall[m][ibin] / (icount[m]*nideal); + vfrac = constant * (rupper*rupper - rlower*rlower); + if (vfrac * normfac != 0.0) + gr = histall[m][ibin] / (vfrac * normfac * icount[m]); else gr = 0.0; - ncoord += gr*nideal; + if (icount[m] != 0) + ncoord += gr * vfrac * normfac; array[ibin][1+2*m] = gr; array[ibin][2+2*m] = ncoord; } diff --git a/src/compute_rdf.h b/src/compute_rdf.h index 794fcc9ee6..9d0956212b 100644 --- a/src/compute_rdf.h +++ b/src/compute_rdf.h @@ -45,6 +45,7 @@ class ComputeRDF : public Compute { int *typecount; int *icount,*jcount; + int *duplicates; class NeighList *list; // half neighbor list }; diff --git a/src/fix_ave_spatial.cpp b/src/fix_ave_spatial.cpp index 521f0956bd..62b7f057af 100644 --- a/src/fix_ave_spatial.cpp +++ b/src/fix_ave_spatial.cpp @@ -888,7 +888,7 @@ void FixAveSpatial::end_of_step() fflush(fp); if (overwrite) { long fileend = ftell(fp); - ftruncate(fileno(fp),fileend); + if (fileend > 0) ftruncate(fileno(fp),fileend); } } } diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index 9c3b9e6753..d502f6a963 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -907,7 +907,7 @@ void FixAveTime::invoke_vector(bigint ntimestep) fflush(fp); if (overwrite) { long fileend = ftell(fp); - ftruncate(fileno(fp),fileend); + if (fileend > 0) ftruncate(fileno(fp),fileend); } } } diff --git a/src/info.cpp b/src/info.cpp index 7b6117acb8..c34eb401cd 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -104,6 +104,7 @@ void Info::command(int narg, char **arg) ++idx; } else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0) && (strncmp(arg[idx+1],"screen",3) == 0)) { + if ((out != screen) && (out != logfile)) fclose(out); out = screen; idx += 2; } else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0) @@ -547,7 +548,7 @@ bool Info::is_available(const char *category, const char *name) int match = 0; if (strcmp(category,"command") == 0) { - if (input->command_map->find(name) != input->command_map->end()); + if (input->command_map->find(name) != input->command_map->end()) match = 1; } else if (strcmp(category,"compute") == 0) { diff --git a/src/pair.cpp b/src/pair.cpp index a87b674608..665615edd0 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -119,6 +119,7 @@ Pair::~Pair() { num_tally_compute = 0; memory->sfree((void *)list_tally_compute); + list_tally_compute = NULL; if (copymode) return; diff --git a/src/velocity.cpp b/src/velocity.cpp index 8b83b007b5..37785af5c3 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -364,7 +364,8 @@ void Velocity::create(double t_desired, int seed) // no-bias compute calculates temp only for new thermal velocities double t; - if (bias_flag == 0) t = temperature->compute_scalar(); + if ((bias_flag == 0) || (temperature_nobias == NULL)) + t = temperature->compute_scalar(); else t = temperature_nobias->compute_scalar(); rescale(t,t_desired);