ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12231 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2014-07-30 14:59:20 +00:00
parent f6c15d997d
commit fe1ade5fad
43 changed files with 379 additions and 156 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/KSPACE/fix_tune_kspace.cpp
View File

@ -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<double>::epsilon()*1.0e-3;
double d=0.0,etemp;

23
src/PERI/pair_peri_eps.cpp
View File

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

51
src/USER-LB/fix_lb_fluid.cpp
View File

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

14
src/USER-OMP/dihedral_quadratic_omp.cpp
View File

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

12
src/USER-OMP/pair_buck_long_coul_long_omp.cpp
View File

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

1
src/USER-OMP/pppm_disp_omp.cpp
View File

@ -59,7 +59,6 @@ void PPPMDispOMP::allocate()
PPPMDisp::allocate();
#if defined(_OPENMP)
const int nthreads = comm->nthreads;
#pragma omp parallel default(none)
#endif
{

2
src/USER-REAXC/fix_qeq_reax.cpp
View File

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

3
src/USER-REAXC/reaxc_hydrogen_bonds.cpp
View File

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

4
src/angle_hybrid.cpp
View File

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

2
src/atom_vec_hybrid.cpp
View File

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

177
src/dump_custom.cpp
View File

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

8
src/dump_custom.h
View File

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

49
src/fix_store_state.cpp
View File

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

26
src/kspace.cpp
View File

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

2
src/lattice.h
View File

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

2
src/library.h
View File

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

2
src/math_complex.h
View File

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

2
src/math_const.h
View File

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

72
src/math_extra.h
View File

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

2
src/math_vector.h
View File

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

2
src/memory.h
View File

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

2
src/min.h
View File

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

2
src/min_cg.h
View File

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

2
src/min_fire.h
View File

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

2
src/min_hftn.h
View File

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

2
src/min_linesearch.h
View File

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

2
src/min_quickmin.h
View File

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

2
src/min_sd.h
View File

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

2
src/minimize.h
View File

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

2
src/modify.h
View File

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

12
src/molecule.cpp
View File

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

2
src/neigh_bond.h
View File

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

2
src/neigh_derive.h
View File

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

2
src/neigh_full.h
View File

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

2
src/neigh_gran.h
View File

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

2
src/neigh_half_bin.h
View File

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

2
src/neigh_half_multi.h
View File

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

2
src/neigh_half_nsq.h
View File

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

2
src/neigh_list.h
View File

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

2
src/neigh_request.h
View File

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

2
src/neigh_respa.h
View File

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

2
src/output.h
View File

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

24
src/pair.cpp
View File

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