From fe1ade5fad2b7e7ca0a8b38627ee2c67bd77435c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 30 Jul 2014 14:59:20 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12231 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/KSPACE/fix_tune_kspace.cpp | 3 +- src/PERI/pair_peri_eps.cpp | 23 +-- src/USER-LB/fix_lb_fluid.cpp | 51 ++--- src/USER-OMP/dihedral_quadratic_omp.cpp | 14 +- src/USER-OMP/pair_buck_long_coul_long_omp.cpp | 12 +- src/USER-OMP/pppm_disp_omp.cpp | 1 - src/USER-REAXC/fix_qeq_reax.cpp | 2 +- src/USER-REAXC/reaxc_hydrogen_bonds.cpp | 3 +- src/angle_hybrid.cpp | 4 + src/atom_vec_hybrid.cpp | 2 + src/dump_custom.cpp | 177 +++++++++++++++++- src/dump_custom.h | 8 +- src/fix_store_state.cpp | 49 ++++- src/kspace.cpp | 26 ++- src/lattice.h | 2 +- src/library.h | 2 +- src/math_complex.h | 2 +- src/math_const.h | 2 + src/math_extra.h | 72 +++---- src/math_vector.h | 2 +- src/memory.h | 2 +- src/min.h | 2 +- src/min_cg.h | 2 +- src/min_fire.h | 2 +- src/min_hftn.h | 2 +- src/min_linesearch.h | 2 +- src/min_quickmin.h | 2 +- src/min_sd.h | 2 +- src/minimize.h | 2 +- src/modify.h | 2 +- src/molecule.cpp | 12 +- src/neigh_bond.h | 2 +- src/neigh_derive.h | 2 +- src/neigh_full.h | 2 +- src/neigh_gran.h | 2 +- src/neigh_half_bin.h | 2 +- src/neigh_half_multi.h | 2 +- src/neigh_half_nsq.h | 2 +- src/neigh_list.h | 2 +- src/neigh_request.h | 2 +- src/neigh_respa.h | 2 +- src/output.h | 2 +- src/pair.cpp | 24 ++- 43 files changed, 379 insertions(+), 156 deletions(-) diff --git a/src/KSPACE/fix_tune_kspace.cpp b/src/KSPACE/fix_tune_kspace.cpp index 0ec6be3e67..9abfc9d1b4 100644 --- a/src/KSPACE/fix_tune_kspace.cpp +++ b/src/KSPACE/fix_tune_kspace.cpp @@ -393,7 +393,7 @@ void FixTuneKspace::adjust_rcut(double time) void FixTuneKspace::mnbrak() { const double GLIMIT = 100.0, TINY = 1.0e-20; - double temp,r,q; + double r,q; r = (bx_brent - ax_brent)*(fb_brent - fc_brent); q = (bx_brent - cx_brent)*(fb_brent - fa_brent); dx_brent = bx_brent - ((bx_brent - cx_brent)*q - (bx_brent - ax_brent)*r)/ @@ -475,7 +475,6 @@ void FixTuneKspace::brent0() void FixTuneKspace::brent1() { - const int ITMAX=100; const double CGOLD=0.3819660; const double ZEPS=numeric_limits::epsilon()*1.0e-3; double d=0.0,etemp; diff --git a/src/PERI/pair_peri_eps.cpp b/src/PERI/pair_peri_eps.cpp index e5b7efbe13..9e06ff9f62 100644 --- a/src/PERI/pair_peri_eps.cpp +++ b/src/PERI/pair_peri_eps.cpp @@ -219,12 +219,10 @@ void PairPeriEPS::compute(int eflag, int vflag) } // ******** temp array to store Plastic extension *********** /// - double deviatorPlasticExtTemp[nlocal][maxpartner]; - for (int ii = 0; ii < nlocal; ii++) { - for (int kk = 0; kk < maxpartner; kk++) { - deviatorPlasticExtTemp[ii][kk] = 0.0; - } - } + // create on heap to reduce stack use and to allow for faster zeroing + double **deviatorPlasticExtTemp; + memory->create(deviatorPlasticExtTemp,nlocal,maxpartner,"pair:plastext"); + memset(&(deviatorPlasticExtTemp[0][0]),0,sizeof(double)*nlocal*maxpartner); // ******** temp array to store Plastic extension *********** /// @@ -402,15 +400,12 @@ void PairPeriEPS::compute(int eflag, int vflag) // store new s0 - for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; + memcpy(s0,s0_new,sizeof(double)*nlocal); - for (i = 0; i < nlocal; i++) { - jnum = npartner[i]; - for (jj = 0; jj < jnum; jj++) { - double temp_data = deviatorPlasticExtTemp[i][jj]; - deviatorPlasticextension[i][jj] = temp_data; - } - } + memcpy(&(deviatorPlasticextension[0][0]), + &(deviatorPlasticExtTemp[0][0]), + sizeof(double)*nlocal*maxpartner); + memory->destroy(deviatorPlasticExtTemp); } /* ---------------------------------------------------------------------- diff --git a/src/USER-LB/fix_lb_fluid.cpp b/src/USER-LB/fix_lb_fluid.cpp index b76f5cc620..4706d576e4 100755 --- a/src/USER-LB/fix_lb_fluid.cpp +++ b/src/USER-LB/fix_lb_fluid.cpp @@ -39,7 +39,6 @@ using namespace LAMMPS_NS; using namespace FixConst; static const double kappa_lb=0.0; -static const double sqrt2=1.41421356237310; FixLbFluid::FixLbFluid(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) @@ -408,69 +407,69 @@ a z wall velocity without implementing fixed BCs in z"); // Create the MPI datatypes used to pass portions of arrays: // datatypes to pass the f and feq arrays. //-------------------------------------------------------------------------- - MPI_Aint sizeofdouble; - MPI_Type_extent(MPI_DOUBLE,&sizeofdouble); + MPI_Aint lb, sizeofdouble; + MPI_Type_get_extent(MPI_DOUBLE,&lb,&sizeofdouble); MPI_Type_vector(subNbz-2,numvel,numvel,MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNby-2,1,numvel*subNbz*sizeofdouble,oneslice,&passxf); + MPI_Type_create_hvector(subNby-2,1,numvel*subNbz*sizeofdouble,oneslice,&passxf); MPI_Type_commit(&passxf); - MPI_Type_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passyf); + MPI_Type_create_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passyf); MPI_Type_commit(&passyf); MPI_Type_free(&oneslice); MPI_Type_vector(subNby,numvel,numvel*subNbz,MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passzf); + MPI_Type_create_hvector(subNbx,1,numvel*subNbz*subNby*sizeofdouble,oneslice,&passzf); MPI_Type_commit(&passzf); // datatypes to pass the u array, and the Ff array. MPI_Type_free(&oneslice); MPI_Type_vector(subNbz+3,3,3,MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxu); + MPI_Type_create_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxu); MPI_Type_commit(&passxu); - MPI_Type_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyu); + MPI_Type_create_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyu); MPI_Type_commit(&passyu); MPI_Type_free(&oneslice); MPI_Type_vector(subNby+3,3,3*(subNbz+3),MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzu); + MPI_Type_create_hvector(subNbx+3,1,3*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzu); MPI_Type_commit(&passzu); // datatypes to pass the density array. MPI_Type_free(&oneslice); MPI_Type_vector(subNbz+3,1,1,MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNby+3,1,1*(subNbz+3)*sizeofdouble,oneslice,&passxrho); + MPI_Type_create_hvector(subNby+3,1,1*(subNbz+3)*sizeofdouble,oneslice,&passxrho); MPI_Type_commit(&passxrho); - MPI_Type_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyrho); + MPI_Type_create_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passyrho); MPI_Type_commit(&passyrho); MPI_Type_free(&oneslice); MPI_Type_vector(subNby+3,1,1*(subNbz+3),MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzrho); + MPI_Type_create_hvector(subNbx+3,1,1*(subNbz+3)*(subNby+3)*sizeofdouble,oneslice,&passzrho); MPI_Type_commit(&passzrho); // datatypes to receive a portion of the Ff array. MPI_Type_free(&oneslice); MPI_Type_vector(subNbz+3,3,3,MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxtemp); + MPI_Type_create_hvector(subNby+3,1,3*(subNbz+3)*sizeofdouble,oneslice,&passxtemp); MPI_Type_commit(&passxtemp); - MPI_Type_hvector(subNbx+3,1,3*(subNbz+3)*5*sizeofdouble,oneslice,&passytemp); + MPI_Type_create_hvector(subNbx+3,1,3*(subNbz+3)*5*sizeofdouble,oneslice,&passytemp); MPI_Type_commit(&passytemp); MPI_Type_free(&oneslice); MPI_Type_vector(subNby+3,3,3*5,MPI_DOUBLE,&oneslice); MPI_Type_commit(&oneslice); - MPI_Type_hvector(subNbx+3,1,3*5*(subNby+3)*sizeofdouble,oneslice,&passztemp); + MPI_Type_create_hvector(subNbx+3,1,3*5*(subNby+3)*sizeofdouble,oneslice,&passztemp); MPI_Type_commit(&passztemp); MPI_Type_free(&oneslice); @@ -806,7 +805,6 @@ void FixLbFluid::calc_fluidforce(void) MPI_Status statuses[20]; double forceloc[3],force[3]; double torqueloc[3],torque[3]; - int numrequests; //-------------------------------------------------------------------------- // Zero out arrays @@ -1317,7 +1315,7 @@ void FixLbFluid::write_restartfile(void) char *hfile; hfile = new char[32]; - sprintf(hfile,"FluidRestart_%d.dat",update->ntimestep); + sprintf(hfile,"FluidRestart_" BIGINT_FORMAT ".dat",update->ntimestep); MPI_File_open(world,hfile,MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL,&fh); @@ -1965,14 +1963,13 @@ void FixLbFluid::initialize_feq(void) void FixLbFluid::equilibriumdist15(int xstart, int xend, int ystart, int yend, int zstart, int zend) { double rho; - double usq; int i, j, k, l, iup, idwn, jup, jdwn, kup, kdwn; double Fx_w, Fy_w, Fz_w; double total_density(0.0); double drhox, drhoy, drhoz, drhoxx, drhoyy, drhozz; double Pxx, Pyy, Pzz, Pxy, Pxz, Pyz; - double grs, p0,TrP; + double grs, p0; double dPdrho; double S[2][3],std; @@ -2064,11 +2061,6 @@ void FixLbFluid::equilibriumdist15(int xstart, int xend, int ystart, int yend, i Fy_w = Ff[i][j][k][1]; Fz_w = Ff[i][j][k][2]; - // this is Tr(P) - TrP = Pxx+Pyy+Pzz; - usq=u_lb[i][j][k][0]*u_lb[i][j][k][0]+u_lb[i][j][k][1]*u_lb[i][j][k][1]+ - u_lb[i][j][k][2]*u_lb[i][j][k][2]; - etacov[0] = rho; etacov[1] = rho*u_lb[i][j][k][0] + Fx_w*tau + rho*bodyforcex*tau; etacov[2] = rho*u_lb[i][j][k][1] + Fy_w*tau + rho*bodyforcey*tau; @@ -2090,6 +2082,7 @@ void FixLbFluid::equilibriumdist15(int xstart, int xend, int ystart, int yend, i etacov[11] = 0.0; etacov[12] = 0.0; etacov[13] = rho*u_lb[i][j][k][0]*u_lb[i][j][k][1]*u_lb[i][j][k][2]; + const double TrP = Pxx+Pyy+Pzz; etacov[14] = K_0*(rho-TrP); for (l=0; l<15; l++) { @@ -2149,14 +2142,13 @@ void FixLbFluid::equilibriumdist15(int xstart, int xend, int ystart, int yend, i void FixLbFluid::equilibriumdist19(int xstart, int xend, int ystart, int yend, int zstart, int zend) { double rho; - double usq; int i, j, k, l, iup, idwn, jup, jdwn, kup, kdwn; double Fx_w, Fy_w, Fz_w; double total_density(0.0); double drhox, drhoy, drhoz, drhoxx, drhoyy, drhozz; double Pxx, Pyy, Pzz, Pxy, Pxz, Pyz; - double grs, p0,TrP; + double grs, p0; double dPdrho; double S[2][3],std; @@ -2247,11 +2239,6 @@ void FixLbFluid::equilibriumdist19(int xstart, int xend, int ystart, int yend, i Fy_w = Ff[i][j][k][1]; Fz_w = Ff[i][j][k][2]; - // this is Tr(P) - TrP = Pxx+Pyy+Pzz; - usq=u_lb[i][j][k][0]*u_lb[i][j][k][0]+u_lb[i][j][k][1]*u_lb[i][j][k][1]+ - u_lb[i][j][k][2]*u_lb[i][j][k][2]; - etacov[0] = rho; etacov[1] = rho*u_lb[i][j][k][0] + Fx_w*tau + rho*bodyforcex*tau; etacov[2] = rho*u_lb[i][j][k][1] + Fy_w*tau + rho*bodyforcey*tau; @@ -2558,8 +2545,6 @@ void FixLbFluid::update_periodic(int xstart, int xend, int ystart, int yend, int void FixLbFluid::streamout(void) { int i,j,k; - double mass,massloc; - double momentumloc[3],momentum[3]; int istart,jstart,kstart; int iend,jend,kend; int w,iproc; diff --git a/src/USER-OMP/dihedral_quadratic_omp.cpp b/src/USER-OMP/dihedral_quadratic_omp.cpp index 067b7bf4e2..2fc72e7317 100644 --- a/src/USER-OMP/dihedral_quadratic_omp.cpp +++ b/src/USER-OMP/dihedral_quadratic_omp.cpp @@ -91,7 +91,7 @@ void DihedralQuadraticOMP::eval(int nfrom, int nto, ThrData * const thr) double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; double c2mag,sc1,sc2,s1,s12,c,p,pd,a,a11,a22; double a33,a12,a13,a23,sx2,sy2,sz2; - double s2,cx,cy,cz,cmag,dx,phi,si,sin2; + double s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; edihedral = 0.0; @@ -212,19 +212,17 @@ void DihedralQuadraticOMP::eval(int nfrom, int nto, ThrData * const thr) // pd = dp/dc phi = acos(c); - if (dx < 0.0) phi *= -1.0; + if (dx > 0.0) phi *= -1.0; si = sin(phi); + if (fabs(si) < SMALLER) si = SMALLER; + siinv = 1.0/si; double dphi = phi-phi0[type]; p = k[type]*dphi; - if (fabs(si) < SMALLER) { - pd = - 2.0 * k[type]; - } else { - pd = - 2.0 * p / si; - } + pd = - 2.0 * p * siinv; p = p * dphi; - if (EFLAG) edihedral = p; + if (EFLAG) edihedral = p; a = pd; c = c * a; diff --git a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp index ce5b7ac3ef..5f5570543b 100644 --- a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp +++ b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairBuckLongCoulLongOMP::PairBuckLongCoulLongOMP(LAMMPS *lmp) : +PairBuckLongCoulLongOMP::PairBuckLongCoulLongOMP(LAMMPS *lmp) : PairBuckLongCoulLong(lmp), ThrOMP(lmp, THR_PAIR) { suffix_flag |= Suffix::OMP; @@ -648,7 +648,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr) // loop over neighbors of my atoms int i, ii, j; - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; + int *jneigh, *jneighn, typei, typej, ni; double qi = 0.0, qri = 0.0, *cutsqi, *cut_bucksqi, *buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti; double r, rsq, r2inv, force_coul, force_buck; @@ -815,7 +815,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t const double cut_out_off_sq = cut_out_off*cut_out_off; - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; + int *jneigh, *jneighn, typei, typej, ni; const int order1 = (ewald_order|(ewald_off^-1))&(1<<1); int i, j, ii; double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi; @@ -1002,7 +1002,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t double *f0 = f[0], *fi = f0; int *ilist = listouter->ilist; - + int i, j, ii; int *jneigh, *jneighn, typei, typej, ni, respa_flag; double qi = 0.0, qri = 0.0; @@ -1011,10 +1011,10 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; double respa_buck = 0.0, respa_coul = 0.0, frespa = 0.0; vector xi, d; - + const double cut_in_off = cut_respa[2]; const double cut_in_on = cut_respa[3]; - + const double cut_in_diff = cut_in_on - cut_in_off; const double cut_in_off_sq = cut_in_off*cut_in_off; const double cut_in_on_sq = cut_in_on*cut_in_on; diff --git a/src/USER-OMP/pppm_disp_omp.cpp b/src/USER-OMP/pppm_disp_omp.cpp index a2ac4f1278..21ba45e91d 100644 --- a/src/USER-OMP/pppm_disp_omp.cpp +++ b/src/USER-OMP/pppm_disp_omp.cpp @@ -59,7 +59,6 @@ void PPPMDispOMP::allocate() PPPMDisp::allocate(); #if defined(_OPENMP) - const int nthreads = comm->nthreads; #pragma omp parallel default(none) #endif { diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index 1af6705ad2..c49f5d98e1 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -48,7 +48,7 @@ using namespace FixConst; //#define LOOSE_ZONE 0.7 #define SQR(x) ((x)*(x)) #define CUBE(x) ((x)*(x)*(x)) -#define MIN_NBRS 100327 +#define MIN_NBRS 100 static const char cite_fix_qeq_reax[] = "fix qeq/reax command:\n\n" diff --git a/src/USER-REAXC/reaxc_hydrogen_bonds.cpp b/src/USER-REAXC/reaxc_hydrogen_bonds.cpp index 640c398fe3..7cbd0bb0a9 100644 --- a/src/USER-REAXC/reaxc_hydrogen_bonds.cpp +++ b/src/USER-REAXC/reaxc_hydrogen_bonds.cpp @@ -42,7 +42,7 @@ void Hydrogen_Bonds( reax_system *system, control_params *control, int itr, top; int num_hb_intrs = 0; ivec rel_jk; - real r_ij, r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; + real r_jk, theta, cos_theta, sin_xhz4, cos_xhz1, sin_theta2; real e_hb, exp_hb2, exp_hb3, CEhb1, CEhb2, CEhb3; rvec dcos_theta_di, dcos_theta_dj, dcos_theta_dk; rvec dvec_jk, force, ext_press; @@ -102,7 +102,6 @@ void Hydrogen_Bonds( reax_system *system, control_params *control, bo_ij = &(pbond_ij->bo_data); type_i = system->my_atoms[i].type; if (type_i < 0) continue; - r_ij = pbond_ij->d; hbp = &(system->reax_param.hbp[ type_i ][ type_j ][ type_k ]); ++num_hb_intrs; diff --git a/src/angle_hybrid.cpp b/src/angle_hybrid.cpp index 784f212653..6f1cceba00 100644 --- a/src/angle_hybrid.cpp +++ b/src/angle_hybrid.cpp @@ -256,6 +256,10 @@ void AngleHybrid::coeff(int narg, char **arg) if (m == nstyles) { if (strcmp(arg[1],"none") == 0) none = 1; else if (strcmp(arg[1],"skip") == 0) none = skip = 1; + else if (strcmp(arg[1],"ba") == 0) + error->all(FLERR,"BondAngle coeff for hybrid angle has invalid format"); + else if (strcmp(arg[1],"bb") == 0) + error->all(FLERR,"BondBond coeff for hybrid angle has invalid format"); else error->all(FLERR,"Angle coeff for hybrid has invalid style"); } diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index 0d82dd327b..9514145d4c 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -109,6 +109,8 @@ void AtomVecHybrid::process_args(int narg, char **arg) dipole_type = MAX(dipole_type,styles[k]->dipole_type); forceclearflag = MAX(forceclearflag,styles[k]->forceclearflag); + if (styles[k]->molecular == 2) onemols = styles[k]->onemols; + comm_x_only = MIN(comm_x_only,styles[k]->comm_x_only); comm_f_only = MIN(comm_f_only,styles[k]->comm_f_only); size_forward += styles[k]->size_forward - 3; diff --git a/src/dump_custom.cpp b/src/dump_custom.cpp index 3cef1beff9..56a38e2402 100644 --- a/src/dump_custom.cpp +++ b/src/dump_custom.cpp @@ -42,7 +42,7 @@ enum{ID,MOL,PROC,PROCP1,TYPE,ELEMENT,MASS, Q,MUX,MUY,MUZ,MU,RADIUS,DIAMETER, OMEGAX,OMEGAY,OMEGAZ,ANGMOMX,ANGMOMY,ANGMOMZ, TQX,TQY,TQZ, - COMPUTE,FIX,VARIABLE}; + COMPUTE,FIX,VARIABLE,INAME,DNAME}; enum{LT,LE,GT,GE,EQ,NEQ}; enum{INT,DOUBLE,STRING,BIGINT}; // same as in DumpCFG @@ -94,6 +94,10 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) : variable = NULL; vbuf = NULL; + ncustom = 0; + id_custom = NULL; + flag_custom = NULL; + // process attributes // ioptional = start of additional optional args // only dump image and dump movie styles process optional args @@ -172,6 +176,10 @@ DumpCustom::~DumpCustom() for (int i = 0; i < nvariable; i++) memory->destroy(vbuf[i]); delete [] vbuf; + for (int i = 0; i < ncustom; i++) delete [] id_custom[i]; + memory->sfree(id_custom); + delete [] flag_custom; + memory->destroy(choose); memory->destroy(dchoose); memory->destroy(clist); @@ -269,6 +277,13 @@ void DumpCustom::init_style() variable[i] = ivariable; } + int icustom; + for (int i = 0; i < ncustom; i++) { + icustom = atom->find_custom(id_custom[i],flag_custom[i]); + if (icustom < 0) + error->all(FLERR,"Could not find custom per-atom property ID"); + } + // set index and check validity of region if (iregion >= 0) { @@ -836,6 +851,24 @@ int DumpCustom::count() i = nfield + ithresh; ptr = vbuf[field2index[i]]; nstride = 1; + + } else if (thresh_array[ithresh] == DNAME) { + int iwhich,tmp; + i = nfield + ithresh; + iwhich = atom->find_custom(id_custom[field2index[i]],tmp); + ptr = atom->dvector[iwhich]; + nstride = 1; + + } else if (thresh_array[ithresh] == INAME) { + int iwhich,tmp; + i = nfield + ithresh; + iwhich = atom->find_custom(id_custom[field2index[i]],tmp); + + int *ivector = atom->ivector[iwhich]; + for (i = 0; i < nlocal; i++) + dchoose[i] = ivector[i]; + ptr = dchoose; + nstride = 1; } // unselect atoms that don't meet threshhold criterion @@ -1248,6 +1281,50 @@ int DumpCustom::parse_fields(int narg, char **arg) field2index[i] = add_variable(suffix); delete [] suffix; + // custom per-atom floating point value = d_ID + + } else if (strncmp(arg[iarg],"d_",2) == 0) { + pack_choice[i] = &DumpCustom::pack_custom; + vtype[i] = DOUBLE; + + int n = strlen(arg[iarg]); + char *suffix = new char[n]; + strcpy(suffix,&arg[iarg][2]); + argindex[i] = 0; + + int tmp = -1; + n = atom->find_custom(suffix,tmp); + if (n < 0) + error->all(FLERR,"Could not find custom per-atom property ID"); + + if (tmp != 1) + error->all(FLERR,"Custom per-atom property ID is not floating point"); + + field2index[i] = add_custom(suffix,1); + delete [] suffix; + + // custom per-atom integer value = i_ID + + } else if (strncmp(arg[iarg],"i_",2) == 0) { + pack_choice[i] = &DumpCustom::pack_custom; + vtype[i] = INT; + + int n = strlen(arg[iarg]); + char *suffix = new char[n]; + strcpy(suffix,&arg[iarg][2]); + argindex[i] = 0; + + int tmp = -1; + n = atom->find_custom(suffix,tmp); + if (n < 0) + error->all(FLERR,"Could not find custom per-atom property ID"); + + if (tmp != 0) + error->all(FLERR,"Custom per-atom property ID is not integer"); + + field2index[i] = add_custom(suffix,0); + delete [] suffix; + } else return iarg; } @@ -1333,6 +1410,34 @@ int DumpCustom::add_variable(char *id) return nvariable-1; } +/* ---------------------------------------------------------------------- + add custom atom property to list used by dump + return index of where this property is in list + if already in list, do not add, just return index, else add to list +------------------------------------------------------------------------- */ + +int DumpCustom::add_custom(char *id, int flag) +{ + int icustom; + for (icustom = 0; icustom < ncustom; icustom++) + if ((strcmp(id,id_custom[icustom]) == 0) + && (flag == flag_custom[icustom])) break; + if (icustom < ncustom) return icustom; + + id_custom = (char **) + memory->srealloc(id_custom,(ncustom+1)*sizeof(char *),"dump:id_custom"); + flag_custom = (int *) + memory->srealloc(flag_custom,(ncustom+1)*sizeof(int),"dump:flag_custom"); + + int n = strlen(id) + 1; + id_custom[ncustom] = new char[n]; + strcpy(id_custom[ncustom],id); + flag_custom[ncustom] = flag; + + ncustom++; + return ncustom-1; +} + /* ---------------------------------------------------------------------- */ int DumpCustom::modify_param(int narg, char **arg) @@ -1574,6 +1679,48 @@ int DumpCustom::modify_param(int narg, char **arg) field2index[nfield+nthresh] = add_variable(suffix); delete [] suffix; + // custom per atom floating point value = d_ID + // must grow field2index and argindex arrays, since access is beyond nfield + + } else if (strncmp(arg[1],"d_",2) == 0) { + thresh_array[nthresh] = DNAME; + memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); + memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); + int n = strlen(arg[1]); + char *suffix = new char[n]; + strcpy(suffix,&arg[1][2]); + argindex[nfield+nthresh] = 0; + + int tmp = -1; + n = atom->find_custom(suffix,tmp); + if ((n < 0) || (tmp != 1)) + error->all(FLERR,"Could not find dump modify " + "custom atom floating point property ID"); + + field2index[nfield+nthresh] = add_custom(suffix,1); + delete [] suffix; + + // custom per atom integer value = i_ID + // must grow field2index and argindex arrays, since access is beyond nfield + + } else if (strncmp(arg[1],"i_",2) == 0) { + thresh_array[nthresh] = INAME; + memory->grow(field2index,nfield+nthresh+1,"dump:field2index"); + memory->grow(argindex,nfield+nthresh+1,"dump:argindex"); + int n = strlen(arg[1]); + char *suffix = new char[n]; + strcpy(suffix,&arg[1][2]); + argindex[nfield+nthresh] = 0; + + int tmp = -1; + n = atom->find_custom(suffix,tmp); + if ((n < 0) || (tmp != 0)) + error->all(FLERR,"Could not find dump modify " + "custom atom integer property ID"); + + field2index[nfield+nthresh] = add_custom(suffix,0); + delete [] suffix; + } else error->all(FLERR,"Invalid dump_modify threshhold operator"); // set operation type of threshhold @@ -1669,6 +1816,34 @@ void DumpCustom::pack_variable(int n) } } +/* ---------------------------------------------------------------------- */ + +void DumpCustom::pack_custom(int n) +{ + + int index = field2index[n]; + + if (flag_custom[index] == 0) { // integer + int iwhich,tmp; + iwhich = atom->find_custom(id_custom[index],tmp); + + int *ivector = atom->ivector[iwhich]; + for (int i = 0; i < nchoose; i++) { + buf[n] = ivector[clist[i]]; + n += size_one; + } + } else if (flag_custom[index] == 1) { // double + int iwhich,tmp; + iwhich = atom->find_custom(id_custom[index],tmp); + + double *dvector = atom->dvector[iwhich]; + for (int i = 0; i < nchoose; i++) { + buf[n] = dvector[clist[i]]; + n += size_one; + } + } +} + /* ---------------------------------------------------------------------- one method for every attribute dump custom can output the atom property is packed into buf starting at n with stride size_one diff --git a/src/dump_custom.h b/src/dump_custom.h index 1fa820ab50..ea3c5ecff8 100644 --- a/src/dump_custom.h +++ b/src/dump_custom.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -69,6 +69,10 @@ class DumpCustom : public Dump { int *variable; // list of indices for the Variables double **vbuf; // local storage for variable evaluation + int ncustom; // # of custom atom properties + char **id_custom; // their names + int *flag_custom; // their data type + int ntypes; // # of atom types char **typenames; // array of element names for each type @@ -86,6 +90,7 @@ class DumpCustom : public Dump { int add_compute(char *); int add_fix(char *); int add_variable(char *); + int add_custom(char *, int); virtual int modify_param(int, char **); typedef void (DumpCustom::*FnPtrHeader)(bigint); @@ -114,6 +119,7 @@ class DumpCustom : public Dump { void pack_compute(int); void pack_fix(int); void pack_variable(int); + void pack_custom(int); void pack_id(int); void pack_molecule(int); diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp index e70eada845..d2ecdbcc5f 100644 --- a/src/fix_store_state.cpp +++ b/src/fix_store_state.cpp @@ -30,7 +30,7 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{KEYWORD,COMPUTE,FIX,VARIABLE}; +enum{KEYWORD,COMPUTE,FIX,VARIABLE,DNAME,INAME}; #define INVOKED_PERATOM 8 @@ -198,11 +198,15 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) : pack_choice[nvalues++] = &FixStoreState::pack_tqz; } else if (strncmp(arg[iarg],"c_",2) == 0 || + strncmp(arg[iarg],"d_",2) == 0 || strncmp(arg[iarg],"f_",2) == 0 || + strncmp(arg[iarg],"i_",2) == 0 || strncmp(arg[iarg],"v_",2) == 0) { cfv_any = 1; if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE; + else if (arg[iarg][0] == 'd') which[nvalues] = DNAME; else if (arg[iarg][0] == 'f') which[nvalues] = FIX; + else if (arg[iarg][0] == 'i') which[nvalues] = INAME; else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE; int n = strlen(arg[iarg]); @@ -265,6 +269,20 @@ FixStoreState::FixStoreState(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Fix store/state compute array is accessed out-of-range"); + } else if (which[i] == INAME) { + int icustom,iflag; + icustom = atom->find_custom(ids[i],iflag); + if ((icustom < 0) || (iflag != 0)) + error->all(FLERR, + "Custom integer vector does not exist"); + + } else if (which[i] == DNAME) { + int icustom,iflag; + icustom = atom->find_custom(ids[i],iflag); + if ((icustom < 0) || (iflag != 1)) + error->all(FLERR, + "Custom floating point vector does not exist"); + } else if (which[i] == FIX) { int ifix = modify->find_fix(ids[i]); if (ifix < 0) @@ -369,6 +387,23 @@ void FixStoreState::init() error->all(FLERR,"Compute ID for fix store/state does not exist"); value2index[m] = icompute; + } else if (which[m] == INAME) { + int icustom,iflag; + icustom = atom->find_custom(ids[m],iflag); + if ((icustom < 0) || (iflag != 0)) + error->all(FLERR, + "Custom integer vector for fix store/state does not exist"); + value2index[m] = icustom; + + } else if (which[m] == DNAME) { + int icustom,iflag; + icustom = atom->find_custom(ids[m],iflag); + if ((icustom < 0) || (iflag != 1)) + error->all(FLERR, + "Custom floating point vector for fix" + " store/state does not exist"); + value2index[m] = icustom; + } else if (which[m] == FIX) { int ifix = modify->find_fix(ids[m]); if (ifix < 0) @@ -465,6 +500,18 @@ void FixStoreState::end_of_step() if (mask[i] & groupbit) values[i][m] = fix_array[i][jm1]; } + // access custom atom property fields + + } else if (which[m] == INAME) { + int *ivector = atom->ivector[n]; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) values[i][m] = ivector[i]; + + } else if (which[m] == DNAME) { + double *dvector = atom->dvector[n]; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) values[i][m] = dvector[i]; + // evaluate atom-style variable } else if (which[m] == VARIABLE) { diff --git a/src/kspace.cpp b/src/kspace.cpp index ac426f4e81..70e0b541d4 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -173,17 +173,31 @@ void KSpace::pair_check() { if (force->pair == NULL) error->all(FLERR,"KSpace solver requires a pair style"); - if (ewaldflag && force->pair->ewaldflag == 0) + + if (ewaldflag && !force->pair->ewaldflag) error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (pppmflag && force->pair->pppmflag == 0) + if (pppmflag && !force->pair->pppmflag) error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (msmflag && force->pair->msmflag == 0) + if (msmflag && !force->pair->msmflag) error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (dispersionflag && force->pair->dispersionflag == 0) + if (dispersionflag && !force->pair->dispersionflag) error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (tip4pflag && force->pair->tip4pflag == 0) + if (tip4pflag && !force->pair->tip4pflag) error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (dipoleflag && force->pair->dipoleflag == 0) + if (dipoleflag && !force->pair->dipoleflag) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + + if (!ewaldflag && force->pair->ewaldflag) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (!pppmflag && force->pair->pppmflag) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (!msmflag && force->pair->msmflag) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (!dispersionflag && force->pair->dispersionflag) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (!tip4pflag && force->pair->tip4pflag) + error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (!dipoleflag && force->pair->dipoleflag) error->all(FLERR,"KSpace style is incompatible with Pair style"); } diff --git a/src/lattice.h b/src/lattice.h index d1ba972fcd..e76fd0507f 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/library.h b/src/library.h index f78c68a3c5..fc85c12809 100644 --- a/src/library.h +++ b/src/library.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/math_complex.h b/src/math_complex.h index 4f53b45288..e4d889c4e6 100644 --- a/src/math_complex.h +++ b/src/math_complex.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/math_const.h b/src/math_const.h index 89437510e3..c0b81bdf8a 100644 --- a/src/math_const.h +++ b/src/math_const.h @@ -25,6 +25,8 @@ namespace MathConst { static const double MY_PI2 = 1.57079632679489661923; // pi/2 static const double MY_PI4 = 0.78539816339744830962; // pi/4 static const double MY_PIS = 1.77245385090551602729; // sqrt(pi) + static const double MY_SQRT2 = 1.41421356237309504880; // sqrt(2) + static const double MY_CBRT2 = 1.25992104989487316476; // 2*(1/3) } } diff --git a/src/math_extra.h b/src/math_extra.h index 8b0af2a460..00186b7056 100755 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -126,7 +126,7 @@ namespace MathExtra { normalize a vector in place ------------------------------------------------------------------------- */ -void MathExtra::norm3(double *v) +inline void MathExtra::norm3(double *v) { double scale = 1.0/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); v[0] *= scale; @@ -138,7 +138,7 @@ void MathExtra::norm3(double *v) normalize a vector, return in ans ------------------------------------------------------------------------- */ -void MathExtra::normalize3(const double *v, double *ans) +inline void MathExtra::normalize3(const double *v, double *ans) { double scale = 1.0/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); ans[0] = v[0]*scale; @@ -150,7 +150,7 @@ void MathExtra::normalize3(const double *v, double *ans) scale a vector to length ------------------------------------------------------------------------- */ -void MathExtra::snormalize3(const double length, const double *v, double *ans) +inline void MathExtra::snormalize3(const double length, const double *v, double *ans) { double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); ans[0] = v[0]*scale; @@ -162,7 +162,7 @@ void MathExtra::snormalize3(const double length, const double *v, double *ans) negate vector v ------------------------------------------------------------------------- */ -void MathExtra::negate3(double *v) +inline void MathExtra::negate3(double *v) { v[0] = -v[0]; v[1] = -v[1]; @@ -173,7 +173,7 @@ void MathExtra::negate3(double *v) scale vector v by s ------------------------------------------------------------------------- */ -void MathExtra::scale3(double s, double *v) +inline void MathExtra::scale3(double s, double *v) { v[0] *= s; v[1] *= s; @@ -184,7 +184,7 @@ void MathExtra::scale3(double s, double *v) ans = v1 + v2 ------------------------------------------------------------------------- */ -void MathExtra::add3(const double *v1, const double *v2, double *ans) +inline void MathExtra::add3(const double *v1, const double *v2, double *ans) { ans[0] = v1[0] + v2[0]; ans[1] = v1[1] + v2[1]; @@ -195,7 +195,7 @@ void MathExtra::add3(const double *v1, const double *v2, double *ans) ans = v1 - v2 ------------------------------------------------------------------------- */ -void MathExtra::sub3(const double *v1, const double *v2, double *ans) +inline void MathExtra::sub3(const double *v1, const double *v2, double *ans) { ans[0] = v1[0] - v2[0]; ans[1] = v1[1] - v2[1]; @@ -206,7 +206,7 @@ void MathExtra::sub3(const double *v1, const double *v2, double *ans) length of vector v ------------------------------------------------------------------------- */ -double MathExtra::len3(const double *v) +inline double MathExtra::len3(const double *v) { return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); } @@ -215,7 +215,7 @@ double MathExtra::len3(const double *v) squared length of vector v, or dot product of v with itself ------------------------------------------------------------------------- */ -double MathExtra::lensq3(const double *v) +inline double MathExtra::lensq3(const double *v) { return v[0]*v[0] + v[1]*v[1] + v[2]*v[2]; } @@ -224,7 +224,7 @@ double MathExtra::lensq3(const double *v) dot product of 2 vectors ------------------------------------------------------------------------- */ -double MathExtra::dot3(const double *v1, const double *v2) +inline double MathExtra::dot3(const double *v1, const double *v2) { return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]; } @@ -233,7 +233,7 @@ double MathExtra::dot3(const double *v1, const double *v2) cross product of 2 vectors ------------------------------------------------------------------------- */ -void MathExtra::cross3(const double *v1, const double *v2, double *ans) +inline void MathExtra::cross3(const double *v1, const double *v2, double *ans) { ans[0] = v1[1]*v2[2] - v1[2]*v2[1]; ans[1] = v1[2]*v2[0] - v1[0]*v2[2]; @@ -262,7 +262,7 @@ void MathExtra::col2mat(const double *ex, const double *ey, const double *ez, determinant of a matrix ------------------------------------------------------------------------- */ -double MathExtra::det3(const double m[3][3]) +inline double MathExtra::det3(const double m[3][3]) { double ans = m[0][0]*m[1][1]*m[2][2] - m[0][0]*m[1][2]*m[2][1] - m[1][0]*m[0][1]*m[2][2] + m[1][0]*m[0][2]*m[2][1] + @@ -274,7 +274,7 @@ double MathExtra::det3(const double m[3][3]) diagonal matrix times a full matrix ------------------------------------------------------------------------- */ -void MathExtra::diag_times3(const double *d, const double m[3][3], +inline void MathExtra::diag_times3(const double *d, const double m[3][3], double ans[3][3]) { ans[0][0] = d[0]*m[0][0]; @@ -310,7 +310,7 @@ void MathExtra::times3_diag(const double m[3][3], const double *d, add two matrices ------------------------------------------------------------------------- */ -void MathExtra::plus3(const double m[3][3], const double m2[3][3], +inline void MathExtra::plus3(const double m[3][3], const double m2[3][3], double ans[3][3]) { ans[0][0] = m[0][0]+m2[0][0]; @@ -328,7 +328,7 @@ void MathExtra::plus3(const double m[3][3], const double m2[3][3], multiply mat1 times mat2 ------------------------------------------------------------------------- */ -void MathExtra::times3(const double m[3][3], const double m2[3][3], +inline void MathExtra::times3(const double m[3][3], const double m2[3][3], double ans[3][3]) { ans[0][0] = m[0][0]*m2[0][0] + m[0][1]*m2[1][0] + m[0][2]*m2[2][0]; @@ -346,7 +346,7 @@ void MathExtra::times3(const double m[3][3], const double m2[3][3], multiply the transpose of mat1 times mat2 ------------------------------------------------------------------------- */ -void MathExtra::transpose_times3(const double m[3][3], const double m2[3][3], +inline void MathExtra::transpose_times3(const double m[3][3], const double m2[3][3], double ans[3][3]) { ans[0][0] = m[0][0]*m2[0][0] + m[1][0]*m2[1][0] + m[2][0]*m2[2][0]; @@ -364,7 +364,7 @@ void MathExtra::transpose_times3(const double m[3][3], const double m2[3][3], multiply mat1 times transpose of mat2 ------------------------------------------------------------------------- */ -void MathExtra::times3_transpose(const double m[3][3], const double m2[3][3], +inline void MathExtra::times3_transpose(const double m[3][3], const double m2[3][3], double ans[3][3]) { ans[0][0] = m[0][0]*m2[0][0] + m[0][1]*m2[0][1] + m[0][2]*m2[0][2]; @@ -383,7 +383,7 @@ void MathExtra::times3_transpose(const double m[3][3], const double m2[3][3], does NOT checks for singular or badly scaled matrix ------------------------------------------------------------------------- */ -void MathExtra::invert3(const double m[3][3], double ans[3][3]) +inline void MathExtra::invert3(const double m[3][3], double ans[3][3]) { double den = m[0][0]*m[1][1]*m[2][2]-m[0][0]*m[1][2]*m[2][1]; den += -m[1][0]*m[0][1]*m[2][2]+m[1][0]*m[0][2]*m[2][1]; @@ -404,7 +404,7 @@ void MathExtra::invert3(const double m[3][3], double ans[3][3]) matrix times vector ------------------------------------------------------------------------- */ -void MathExtra::matvec(const double m[3][3], const double *v, double *ans) +inline void MathExtra::matvec(const double m[3][3], const double *v, double *ans) { ans[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2]; ans[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2]; @@ -415,7 +415,7 @@ void MathExtra::matvec(const double m[3][3], const double *v, double *ans) matrix times vector ------------------------------------------------------------------------- */ -void MathExtra::matvec(const double *ex, const double *ey, const double *ez, +inline void MathExtra::matvec(const double *ex, const double *ey, const double *ez, const double *v, double *ans) { ans[0] = ex[0]*v[0] + ey[0]*v[1] + ez[0]*v[2]; @@ -427,7 +427,7 @@ void MathExtra::matvec(const double *ex, const double *ey, const double *ez, transposed matrix times vector ------------------------------------------------------------------------- */ -void MathExtra::transpose_matvec(const double m[3][3], const double *v, +inline void MathExtra::transpose_matvec(const double m[3][3], const double *v, double *ans) { ans[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2]; @@ -439,7 +439,7 @@ void MathExtra::transpose_matvec(const double m[3][3], const double *v, transposed matrix times vector ------------------------------------------------------------------------- */ -void MathExtra::transpose_matvec(const double *ex, const double *ey, +inline void MathExtra::transpose_matvec(const double *ex, const double *ey, const double *ez, const double *v, double *ans) { @@ -452,7 +452,7 @@ void MathExtra::transpose_matvec(const double *ex, const double *ey, transposed matrix times diagonal matrix ------------------------------------------------------------------------- */ -void MathExtra::transpose_diag3(const double m[3][3], const double *d, +inline void MathExtra::transpose_diag3(const double m[3][3], const double *d, double ans[3][3]) { ans[0][0] = m[0][0]*d[0]; @@ -470,7 +470,7 @@ void MathExtra::transpose_diag3(const double m[3][3], const double *d, row vector times matrix ------------------------------------------------------------------------- */ -void MathExtra::vecmat(const double *v, const double m[3][3], double *ans) +inline void MathExtra::vecmat(const double *v, const double m[3][3], double *ans) { ans[0] = v[0]*m[0][0] + v[1]*m[1][0] + v[2]*m[2][0]; ans[1] = v[0]*m[0][1] + v[1]*m[1][1] + v[2]*m[2][1]; @@ -493,7 +493,7 @@ inline void MathExtra::scalar_times3(const double f, double m[3][3]) upper-triangular 3x3, stored as 6-vector in Voigt notation ------------------------------------------------------------------------- */ -void MathExtra::multiply_shape_shape(const double *one, const double *two, +inline void MathExtra::multiply_shape_shape(const double *one, const double *two, double *ans) { ans[0] = one[0]*two[0]; @@ -508,7 +508,7 @@ void MathExtra::multiply_shape_shape(const double *one, const double *two, normalize a quaternion ------------------------------------------------------------------------- */ -void MathExtra::qnormalize(double *q) +inline void MathExtra::qnormalize(double *q) { double norm = 1.0 / sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); q[0] *= norm; @@ -522,7 +522,7 @@ void MathExtra::qnormalize(double *q) assume q is of unit length ------------------------------------------------------------------------- */ -void MathExtra::qconjugate(double *q, double *qc) +inline void MathExtra::qconjugate(double *q, double *qc) { qc[0] = q[0]; qc[1] = -q[1]; @@ -534,7 +534,7 @@ void MathExtra::qconjugate(double *q, double *qc) vector-quaternion multiply: c = a*b, where a = (0,a) ------------------------------------------------------------------------- */ -void MathExtra::vecquat(double *a, double *b, double *c) +inline void MathExtra::vecquat(double *a, double *b, double *c) { c[0] = -a[0]*b[1] - a[1]*b[2] - a[2]*b[3]; c[1] = b[0]*a[0] + a[1]*b[3] - a[2]*b[2]; @@ -546,7 +546,7 @@ void MathExtra::vecquat(double *a, double *b, double *c) quaternion-vector multiply: c = a*b, where b = (0,b) ------------------------------------------------------------------------- */ -void MathExtra::quatvec(double *a, double *b, double *c) +inline void MathExtra::quatvec(double *a, double *b, double *c) { c[0] = -a[1]*b[0] - a[2]*b[1] - a[3]*b[2]; c[1] = a[0]*b[0] + a[2]*b[2] - a[3]*b[1]; @@ -558,7 +558,7 @@ void MathExtra::quatvec(double *a, double *b, double *c) quaternion-quaternion multiply: c = a*b ------------------------------------------------------------------------- */ -void MathExtra::quatquat(double *a, double *b, double *c) +inline void MathExtra::quatquat(double *a, double *b, double *c) { c[0] = a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3]; c[1] = a[0]*b[1] + b[0]*a[1] + a[2]*b[3] - a[3]*b[2]; @@ -573,7 +573,7 @@ void MathExtra::quatquat(double *a, double *b, double *c) c is a three component vector ------------------------------------------------------------------------- */ -void MathExtra::invquatvec(double *a, double *b, double *c) +inline void MathExtra::invquatvec(double *a, double *b, double *c) { c[0] = -a[1]*b[0] + a[0]*b[1] + a[3]*b[2] - a[2]*b[3]; c[1] = -a[2]*b[0] - a[3]*b[1] + a[0]*b[2] + a[1]*b[3]; @@ -585,7 +585,7 @@ void MathExtra::invquatvec(double *a, double *b, double *c) v MUST be a unit vector ------------------------------------------------------------------------- */ -void MathExtra::axisangle_to_quat(const double *v, const double angle, +inline void MathExtra::axisangle_to_quat(const double *v, const double angle, double *quat) { double halfa = 0.5*angle; @@ -600,7 +600,7 @@ void MathExtra::axisangle_to_quat(const double *v, const double angle, Apply principal rotation generator about x to rotation matrix m ------------------------------------------------------------------------- */ -void MathExtra::rotation_generator_x(const double m[3][3], double ans[3][3]) +inline void MathExtra::rotation_generator_x(const double m[3][3], double ans[3][3]) { ans[0][0] = 0; ans[0][1] = -m[0][2]; @@ -617,7 +617,7 @@ void MathExtra::rotation_generator_x(const double m[3][3], double ans[3][3]) Apply principal rotation generator about y to rotation matrix m ------------------------------------------------------------------------- */ -void MathExtra::rotation_generator_y(const double m[3][3], double ans[3][3]) +inline void MathExtra::rotation_generator_y(const double m[3][3], double ans[3][3]) { ans[0][0] = m[0][2]; ans[0][1] = 0; @@ -634,7 +634,7 @@ void MathExtra::rotation_generator_y(const double m[3][3], double ans[3][3]) Apply principal rotation generator about z to rotation matrix m ------------------------------------------------------------------------- */ -void MathExtra::rotation_generator_z(const double m[3][3], double ans[3][3]) +inline void MathExtra::rotation_generator_z(const double m[3][3], double ans[3][3]) { ans[0][0] = -m[0][1]; ans[0][1] = m[0][0]; diff --git a/src/math_vector.h b/src/math_vector.h index da65e069bd..000c2fc07f 100644 --- a/src/math_vector.h +++ b/src/math_vector.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/memory.h b/src/memory.h index bcd59ab68f..b14d0590c2 100644 --- a/src/memory.h +++ b/src/memory.h @@ -408,7 +408,7 @@ class Memory : protected Pointers { TYPE ***plane = (TYPE ***) smalloc(nbytes,name); nbytes = ((bigint) sizeof(TYPE ***)) * n1; array = (TYPE ****) smalloc(nbytes,name); - + int i,j,k; bigint m1,m2; bigint n = 0; diff --git a/src/min.h b/src/min.h index fc78fd8735..417cef1c79 100644 --- a/src/min.h +++ b/src/min.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/min_cg.h b/src/min_cg.h index 6ac8901d24..bf199408e3 100644 --- a/src/min_cg.h +++ b/src/min_cg.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/min_fire.h b/src/min_fire.h index 9f1be75300..3cb375163e 100644 --- a/src/min_fire.h +++ b/src/min_fire.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/min_hftn.h b/src/min_hftn.h index 9d9932572b..c96834a50f 100644 --- a/src/min_hftn.h +++ b/src/min_hftn.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/min_linesearch.h b/src/min_linesearch.h index 51e1277506..27e34fe6eb 100644 --- a/src/min_linesearch.h +++ b/src/min_linesearch.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/min_quickmin.h b/src/min_quickmin.h index 443961972a..a4766bbcdc 100644 --- a/src/min_quickmin.h +++ b/src/min_quickmin.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/min_sd.h b/src/min_sd.h index 2104280729..705103b4d2 100644 --- a/src/min_sd.h +++ b/src/min_sd.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/minimize.h b/src/minimize.h index 2f41291563..4c44618948 100644 --- a/src/minimize.h +++ b/src/minimize.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/modify.h b/src/modify.h index 9761af5200..422c77d2fe 100644 --- a/src/modify.h +++ b/src/modify.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/molecule.cpp b/src/molecule.cpp index a7fa72af83..af2dfaf60a 100644 --- a/src/molecule.cpp +++ b/src/molecule.cpp @@ -718,21 +718,21 @@ void Molecule::dihedrals(int flag, char *line) dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; - dihedral_atom4[m][num_dihedral[m]] = atom4; + dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; m = atom3-1; dihedral_type[m][num_dihedral[m]] = itype; dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; - dihedral_atom4[m][num_dihedral[m]] = atom4; + dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; m = atom4-1; dihedral_type[m][num_dihedral[m]] = itype; dihedral_atom1[m][num_dihedral[m]] = atom1; dihedral_atom2[m][num_dihedral[m]] = atom2; dihedral_atom3[m][num_dihedral[m]] = atom3; - dihedral_atom4[m][num_dihedral[m]] = atom4; + dihedral_atom4[m][num_dihedral[m]] = atom4; num_dihedral[m]++; } } else { @@ -803,21 +803,21 @@ void Molecule::impropers(int flag, char *line) improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; - improper_atom4[m][num_improper[m]] = atom4; + improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; m = atom3-1; improper_type[m][num_improper[m]] = itype; improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; - improper_atom4[m][num_improper[m]] = atom4; + improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; m = atom4-1; improper_type[m][num_improper[m]] = itype; improper_atom1[m][num_improper[m]] = atom1; improper_atom2[m][num_improper[m]] = atom2; improper_atom3[m][num_improper[m]] = atom3; - improper_atom4[m][num_improper[m]] = atom4; + improper_atom4[m][num_improper[m]] = atom4; num_improper[m]++; } } else { diff --git a/src/neigh_bond.h b/src/neigh_bond.h index c9a1e9317a..894c4aebae 100644 --- a/src/neigh_bond.h +++ b/src/neigh_bond.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_derive.h b/src/neigh_derive.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_derive.h +++ b/src/neigh_derive.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_full.h b/src/neigh_full.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_full.h +++ b/src/neigh_full.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_gran.h b/src/neigh_gran.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_gran.h +++ b/src/neigh_gran.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_half_bin.h b/src/neigh_half_bin.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_half_bin.h +++ b/src/neigh_half_bin.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_half_multi.h b/src/neigh_half_multi.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_half_multi.h +++ b/src/neigh_half_multi.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_half_nsq.h b/src/neigh_half_nsq.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_half_nsq.h +++ b/src/neigh_half_nsq.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_list.h b/src/neigh_list.h index f7fecce9b8..2374fde652 100644 --- a/src/neigh_list.h +++ b/src/neigh_list.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_request.h b/src/neigh_request.h index f41a4017b9..769d5354bf 100644 --- a/src/neigh_request.h +++ b/src/neigh_request.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/neigh_respa.h b/src/neigh_respa.h index 5ad76c98af..1538f7662a 100644 --- a/src/neigh_respa.h +++ b/src/neigh_respa.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/output.h b/src/output.h index 9274fdd7ee..cda9e800ee 100644 --- a/src/output.h +++ b/src/output.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov diff --git a/src/pair.cpp b/src/pair.cpp index 44220e276a..fb9add2078 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -42,7 +42,7 @@ using namespace LAMMPS_NS; #define EWALD_F 1.12837917 -enum{RLINEAR,RSQ,BMP}; +enum{NONE,RLINEAR,RSQ,BMP}; /* ---------------------------------------------------------------------- */ @@ -638,15 +638,14 @@ void Pair::free_disp_tables() double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2) { - double value; if (mix_flag == GEOMETRIC) - value = sqrt(eps1*eps2); + return sqrt(eps1*eps2); else if (mix_flag == ARITHMETIC) - value = sqrt(eps1*eps2); + return sqrt(eps1*eps2); else if (mix_flag == SIXTHPOWER) - value = 2.0 * sqrt(eps1*eps2) * - pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0)); - return value; + return (2.0 * sqrt(eps1*eps2) * + pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0))); + else return 0.0; } /* ---------------------------------------------------------------------- @@ -655,14 +654,13 @@ double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2) double Pair::mix_distance(double sig1, double sig2) { - double value; if (mix_flag == GEOMETRIC) - value = sqrt(sig1*sig2); + return sqrt(sig1*sig2); else if (mix_flag == ARITHMETIC) - value = 0.5 * (sig1+sig2); + return (0.5 * (sig1+sig2)); else if (mix_flag == SIXTHPOWER) - value = pow((0.5 * (pow(sig1,6.0) + pow(sig2,6.0))),1.0/6.0); - return value; + return pow((0.5 * (pow(sig1,6.0) + pow(sig2,6.0))),1.0/6.0); + else return 0.0; } /* ---------------------------------------------------------------------- */ @@ -1470,7 +1468,7 @@ void Pair::write_file(int narg, char **arg) int n = force->inumeric(FLERR,arg[2]); - int style; + int style = NONE; if (strcmp(arg[3],"r") == 0) style = RLINEAR; else if (strcmp(arg[3],"rsq") == 0) style = RSQ; else if (strcmp(arg[3],"bitmap") == 0) style = BMP;