diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 925e8235b2..79f49607fd 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -43,6 +43,7 @@ Ewald::Ewald(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg != 1) error->all(FLERR,"Illegal kspace_style ewald command"); + ewaldflag = 1; group_group_enable = 1; group_allocate_flag = 0; @@ -104,8 +105,8 @@ void Ewald::init() scale = 1.0; - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); + pair_check(); + int itmp; double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp); if (p_cutoff == NULL) diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp index d5bb008194..8673102b30 100644 --- a/src/KSPACE/ewald_disp.cpp +++ b/src/KSPACE/ewald_disp.cpp @@ -65,6 +65,8 @@ EwaldDisp::EwaldDisp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) b2 = 0; } +/* ---------------------------------------------------------------------- */ + EwaldDisp::~EwaldDisp() { deallocate(); @@ -73,7 +75,6 @@ EwaldDisp::~EwaldDisp() delete [] B; } - /* --------------------------------------------------------------------- */ void EwaldDisp::init() @@ -88,6 +89,8 @@ void EwaldDisp::init() if (logfile) fprintf(logfile,"EwaldDisp initialization ...\n"); } + if (force->pair == NULL || force->pair->ewaldflag == 0) + error->all(FLERR,"KSpace style is incompatible with Pair style"); if (domain->dimension == 2) // check for errors error->all(FLERR,"Cannot use EwaldDisp with 2d simulation"); if (slabflag == 0 && domain->nonperiodic > 0) @@ -103,6 +106,8 @@ void EwaldDisp::init() //dielectric = force->dielectric; mumurd2e = dielectric = 1.0; + pair_check(); + int tmp; Pair *pair = force->pair; int *ptr = pair ? (int *) pair->extract("ewald_order",tmp) : NULL; @@ -134,7 +139,6 @@ void EwaldDisp::init() nsums += n[k]; } - g_ewald = 0; pair->init(); // so B is defined init_coeffs(); @@ -146,16 +150,12 @@ void EwaldDisp::init() qsum = sum[0].x; qsqsum = sum[0].x2; } - if (function[1]) { - bsbsum = sum[1].x2; - } - if (function[2]) { - bsbsum = sum[2].x2; - } - + if (function[1]) bsbsum = sum[1].x2; + if (function[2]) bsbsum = sum[2].x2; if (qsqsum == 0.0 && bsbsum == 0.0) - error->all(FLERR,"Cannot use Ewald/n solver on system with no charge or LJ particles"); + error->all(FLERR,"Cannot use Ewald/disp solver " + "on system with no charge or LJ particles"); if (fabs(qsum) > SMALL && comm->me == 0) { char str[128]; sprintf(str,"System is not charge neutral, net charge = %g",qsum); @@ -173,31 +173,28 @@ void EwaldDisp::init() b2 = bsbsum; //Are these units right? bigint natoms = atom->natoms; - if (function[0]) { //Coulombic + if (function[0]) { g_ewald = accuracy*sqrt(natoms*(*cutoff)*shape_det(domain->h)) / (2.0*q2); if (g_ewald >= 1.0) error->all(FLERR,"KSpace accuracy too large to estimate G vector"); g_ewald = sqrt(-log(g_ewald)) / *cutoff; } - else if (function[1] || function[2]) { //Only LJ - + else if (function[1] || function[2]) { double *cutoffLJ = pair ? (double *) pair->extract("cut_LJ",tmp) : NULL; - //Try Newton Solver - //Use old method to get guess g_ewald = (1.35 - 0.15*log(accuracy))/ *cutoffLJ; - double g_ewald_new = NewtonSolve(g_ewald,(*cutoffLJ),natoms,shape_det(domain->h),b2); - + double g_ewald_new = + NewtonSolve(g_ewald,(*cutoffLJ),natoms,shape_det(domain->h),b2); if (g_ewald_new > 0.0) g_ewald = g_ewald_new; - else error->warning(FLERR,"Ewald/n Newton solver failed, using old method to estimate g_ewald"); - + else error->warning(FLERR,"Ewald/disp Newton solver failed, " + "using old method to estimate g_ewald"); if (g_ewald >= 1.0) error->all(FLERR,"KSpace accuracy too large to estimate G vector"); } - if (!comm->me) { // output results + if (!comm->me) { if (screen) fprintf(screen, " G vector = %g\n", g_ewald); if (logfile) fprintf(logfile, " G vector = %g\n", g_ewald); } @@ -214,12 +211,13 @@ void EwaldDisp::init() void EwaldDisp::setup() { - volume = shape_det(domain->h)*slab_volfactor; // cell volume - memcpy(unit, domain->h_inv, sizeof(shape)); // wave vector units + volume = shape_det(domain->h)*slab_volfactor; + memcpy(unit, domain->h_inv, sizeof(shape)); shape_scalar_mult(unit, 2.0*MY_PI); unit[2] /= slab_volfactor; //int nbox_old = nbox, nkvec_old = nkvec; + if (accuracy>=1) { nbox = 0; error->all(FLERR,"KSpace accuracy too low"); @@ -255,12 +253,12 @@ void EwaldDisp::setup() gsqmx *= 1.00001; reallocate(); - coefficients(); // compute coeffs + coefficients(); init_coeffs(); init_coeff_sums(); init_self(); - if (!(first_output||comm->me)) { // output on first + if (!(first_output||comm->me)) { first_output = 1; if (screen) fprintf(screen, " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); @@ -277,7 +275,7 @@ double EwaldDisp::rms(int km, double prd, bigint natoms, double q2, double b2) { double value = 0.0; - //Coulombic + // Coulombic double g2 = g_ewald*g_ewald; @@ -285,7 +283,7 @@ double EwaldDisp::rms(int km, double prd, bigint natoms, double q2, double b2) sqrt(1.0/(MY_PI*km*natoms)) * exp(-MY_PI*MY_PI*km*km/(g2*prd*prd)); - //Lennard-Jones + // Lennard-Jones double g7 = g2*g2*g2*g_ewald; @@ -297,13 +295,13 @@ double EwaldDisp::rms(int km, double prd, bigint natoms, double q2, double b2) return value; } -void EwaldDisp::reallocate() // allocate memory +void EwaldDisp::reallocate() { int ix, iy, iz; int nkvec_max = nkvec; vector h; - nkvec = 0; // determine size(kvec) + nkvec = 0; int *kflag = new int[(nbox+1)*(2*nbox+1)*(2*nbox+1)]; int *flag = kflag; @@ -337,8 +335,8 @@ void EwaldDisp::reallocate() // allocate memory nkvec_max = nkvec; } - flag = kflag; // create index and - kvector *k = kvec; // wave vectors + flag = kflag; // create index and + kvector *k = kvec; // wave vectors hvector *hi = hvec; for (ix=0; ix<=nbox; ++ix) for (iy=-nbox; iy<=nbox; ++iy) @@ -397,7 +395,7 @@ void EwaldDisp::deallocate() // free memory } -void EwaldDisp::coefficients() // set up pre-factors +void EwaldDisp::coefficients() { vector h; hvector *hi = hvec, *nh; @@ -442,8 +440,7 @@ void EwaldDisp::coefficients() // set up pre-fact } } - -void EwaldDisp::init_coeffs() // local pair coeffs +void EwaldDisp::init_coeffs() { int tmp; int n = atom->ntypes; @@ -476,10 +473,9 @@ void EwaldDisp::init_coeffs() // local pair coeff } } - -void EwaldDisp::init_coeff_sums() // sums based on atoms +void EwaldDisp::init_coeff_sums() { - if (sums) return; // calculated only once + if (sums) return; // calculated only once sums = 1; Sum sum_local[EWALD_MAX_NSUMS]; @@ -679,7 +675,7 @@ void EwaldDisp::compute_ek() C_ANGLE(z1.y, *(x++)*unit[1]); C_ANGLE(z1.z, *(x++)*unit[2]); } - for (; zzx, zz->x, z1.x); // 3D k-vector C_RMULT(zy->y, zz->y, z1.y); C_CONJ(zx->y, zy->y); C_RMULT(zy->z, zz->z, z1.z); C_CONJ(zx->z, zy->z); @@ -696,7 +692,7 @@ void EwaldDisp::compute_ek() h = hvec; } for (k=kvec; ky) { // based on order in + if (ky!=k->y) { // based on order in if (kx!=k->x) cx = z[kx = k->x].x; // reallocate C_RMULT(zxy, z[ky = k->y].y, cx); } @@ -744,8 +740,8 @@ void EwaldDisp::compute_force() if (atom->torque) t = atom->torque[0]; memcpy(func, function, EWALD_NFUNCS*sizeof(int)); memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector)); // fj = -dE/dr = - for (; fy) { // based on order in + if (ky!=k->y) { // based on order in if (kx!=k->x) zx = z[kx = k->x].x; // reallocate C_RMULT(zxy, z[ky = k->y].y, zx); } @@ -861,7 +857,7 @@ void EwaldDisp::compute_energy() memcpy(func, function, EWALD_NFUNCS*sizeof(int)); memset(sum, 0, EWALD_NFUNCS*sizeof(double)); // reset sums - for (int k=0; kre*cek->re+cek->im*cek->im); ++cek; } if (func[1]) { // geometric 1/r^6 @@ -916,7 +912,7 @@ void EwaldDisp::compute_energy_peratom() mu++; } for (nh = (h = hvec)+nkvec; hy) { // based on order in + if (ky!=k->y) { // based on order in if (kx!=k->x) zx = z[kx = k->x].x; // reallocate C_RMULT(zxy, z[ky = k->y].y, zx); } @@ -981,7 +977,7 @@ void EwaldDisp::compute_virial() memcpy(func, function, EWALD_NFUNCS*sizeof(int)); memset(sum, 0, EWALD_NFUNCS*sizeof(shape)); - for (int k=0; kre*cek->re+cek->im*cek->im; ++cek; sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r; @@ -1050,7 +1046,7 @@ void EwaldDisp::compute_virial_peratom() mu++; } for (nh = (h = hvec)+nkvec; hy) { // based on order in + if (ky!=k->y) { // based on order in if (kx!=k->x) zx = z[kx = k->x].x; // reallocate C_RMULT(zxy, z[ky = k->y].y, zx); } @@ -1180,7 +1176,8 @@ void EwaldDisp::compute_slabcorr() Newton solver used to find g_ewald for LJ systems ------------------------------------------------------------------------- */ -double EwaldDisp::NewtonSolve(double x, double Rc, bigint natoms, double vol, double b2) +double EwaldDisp::NewtonSolve(double x, double Rc, + bigint natoms, double vol, double b2) { double dx,tol; int maxit; @@ -1214,7 +1211,8 @@ double EwaldDisp::f(double x, double Rc, bigint natoms, double vol, double b2) Calculate numerical derivative f'(x) ------------------------------------------------------------------------- */ -double EwaldDisp::derivf(double x, double Rc, bigint natoms, double vol, double b2) +double EwaldDisp::derivf(double x, double Rc, + bigint natoms, double vol, double b2) { double h = 0.000001; //Derivative step-size return (f(x + h,Rc,natoms,vol,b2) - f(x,Rc,natoms,vol,b2)) / h; diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 5a77ce4085..2ce66b673d 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -52,6 +52,8 @@ MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg < 1) error->all(FLERR,"Illegal kspace_style msm command"); + msmflag = 1; + accuracy_relative = atof(arg[0]); nfactors = 1; @@ -80,7 +82,6 @@ MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) levels = 0; differentiation_flag = 1; - } /* ---------------------------------------------------------------------- @@ -154,19 +155,14 @@ void MSM::init() qqrd2e = force->qqrd2e; scale = 1.0; - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); + pair_check(); + int itmp; double *p_cutoff = (double *) force->pair->extract("cut_msm",itmp); if (p_cutoff == NULL) error->all(FLERR,"KSpace style is incompatible with Pair style"); cutoff = *p_cutoff; - if ((strcmp(force->kspace_style,"pppm/tip4p") == 0) || - (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0)) { - error->all(FLERR,"KSpace style is incompatible with Pair style"); - } - // compute qsum & qsqsum and give error if not charge-neutral qsum = qsqsum = 0.0; @@ -184,10 +180,13 @@ void MSM::init() if (qsqsum == 0.0) error->all(FLERR,"Cannot use kspace solver on system with no charge"); + + // not yet sure of the correction needed for non-neutral systems + if (fabs(qsum) > SMALL) { char str[128]; sprintf(str,"System is not charge neutral, net charge = %g",qsum); - error->all(FLERR,str); // Not yet sure of the correction needed for non-neutral systems + error->all(FLERR,str); } // set accuracy (force units) from accuracy_relative or accuracy_absolute diff --git a/src/KSPACE/pair_born_coul_long.cpp b/src/KSPACE/pair_born_coul_long.cpp index e1eb08d20e..498b0099a3 100644 --- a/src/KSPACE/pair_born_coul_long.cpp +++ b/src/KSPACE/pair_born_coul_long.cpp @@ -43,7 +43,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairBornCoulLong::PairBornCoulLong(LAMMPS *lmp) : Pair(lmp) {} +PairBornCoulLong::PairBornCoulLong(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; +} /* ---------------------------------------------------------------------- */ @@ -349,7 +352,7 @@ void PairBornCoulLong::init_style() // insure use of KSpace long-range solver, set g_ewald if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); + error->all(FLERR,"Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; neighbor->request(this); diff --git a/src/KSPACE/pair_born_coul_msm.cpp b/src/KSPACE/pair_born_coul_msm.cpp index 85d960145d..3590265d8b 100644 --- a/src/KSPACE/pair_born_coul_msm.cpp +++ b/src/KSPACE/pair_born_coul_msm.cpp @@ -35,7 +35,11 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp) {} +PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp) +{ + ewaldflag = pppmflag = 0; + msmflag = 1; +} /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_buck_coul_long.cpp b/src/KSPACE/pair_buck_coul_long.cpp index 393cfc1b9b..a29e58aa89 100644 --- a/src/KSPACE/pair_buck_coul_long.cpp +++ b/src/KSPACE/pair_buck_coul_long.cpp @@ -39,7 +39,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairBuckCoulLong::PairBuckCoulLong(LAMMPS *lmp) : Pair(lmp) {} +PairBuckCoulLong::PairBuckCoulLong(LAMMPS *lmp) : Pair(lmp) +{ + ewaldflag = pppmflag = 1; +} /* ---------------------------------------------------------------------- */ @@ -222,7 +225,8 @@ void PairBuckCoulLong::settings(int narg, char **arg) void PairBuckCoulLong::coeff(int narg, char **arg) { - if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 5 || narg > 6) + error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -326,7 +330,7 @@ void PairBuckCoulLong::init_style() // insure use of KSpace long-range solver, set g_ewald if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); + error->all(FLERR,"Pair style requres a KSpace style"); g_ewald = force->kspace->g_ewald; neighbor->request(this); diff --git a/src/KSPACE/pair_buck_coul_msm.cpp b/src/KSPACE/pair_buck_coul_msm.cpp index ffaf111992..22d3e1a64f 100644 --- a/src/KSPACE/pair_buck_coul_msm.cpp +++ b/src/KSPACE/pair_buck_coul_msm.cpp @@ -32,7 +32,11 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp) {} +PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp) +{ + ewaldflag = pppmflag = 0; + msmflag = 1; +} /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_buck_disp_coul_long.cpp b/src/KSPACE/pair_buck_disp_coul_long.cpp index 1c6d0c07e8..15f9023636 100644 --- a/src/KSPACE/pair_buck_disp_coul_long.cpp +++ b/src/KSPACE/pair_buck_disp_coul_long.cpp @@ -48,6 +48,7 @@ using namespace LAMMPS_NS; PairBuckDispCoulLong::PairBuckDispCoulLong(LAMMPS *lmp) : Pair(lmp) { + dispersionflag = ewaldflag = pppmflag = 1; respa_enable = 1; ftable = NULL; } @@ -64,21 +65,24 @@ PairBuckDispCoulLong::PairBuckDispCoulLong(LAMMPS *lmp) : Pair(lmp) #define PAIR_LARGEST "Using largest cut-off for buck/coul long long" #define PAIR_MIX "Geometric mixing assumed for 1/r^6 coefficients" + + void PairBuckDispCoulLong::options(char **arg, int order) { - char *option[] = {"long", "cut", "off", NULL}; + const char *option[] = {"long", "cut", "off", NULL}; int i; if (!*arg) error->all(FLERR,PAIR_ILLEGAL); for (i=0; option[i]&&strcmp(arg[0], option[i]); ++i); switch (i) { default: error->all(FLERR,PAIR_ILLEGAL); - case 0: ewald_order |= 1<me && ewald_order&(1<<6)) error->warning(FLERR,PAIR_MIX); - if (!comm->me && ewald_order==((1<<1)|(1<<6))) error->warning(FLERR,PAIR_LARGEST); + if (!comm->me && ewald_order==((1<<1)|(1<<6))) + error->warning(FLERR,PAIR_LARGEST); if (!*(++arg)) error->all(FLERR,PAIR_MISSING); if (ewald_off&(1<<6)) error->all(FLERR,PAIR_LJ_OFF); if (!((ewald_order^ewald_off)&(1<<1))) error->all(FLERR,PAIR_COUL_CUT); @@ -98,7 +105,7 @@ void PairBuckDispCoulLong::settings(int narg, char **arg) if (narg == 4) cut_coul = force->numeric(*arg); else cut_coul = cut_buck_global; - if (allocated) { // reset explicit cuts + if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) @@ -189,7 +196,8 @@ void *PairBuckDispCoulLong::extract(const char *id, int &dim) void PairBuckDispCoulLong::coeff(int narg, char **arg) { - if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients"); + if (narg < 5 || narg > 6) + error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; @@ -224,7 +232,6 @@ void PairBuckDispCoulLong::coeff(int narg, char **arg) void PairBuckDispCoulLong::init_style() { - // require an atom style with charge defined if (!atom->q_flag && (ewald_order&(1<<1))) @@ -276,18 +283,12 @@ void PairBuckDispCoulLong::init_style() cut_respa = ((Respa *) update->integrate)->cutoff; else cut_respa = NULL; - // ensure use of KSpace long-range solver, set g_ewald + // ensure use of KSpace long-range solver, set two g_ewalds - if (ewald_order&(1<<1)) { // r^-1 kspace - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - g_ewald = force->kspace->g_ewald; - } - if (ewald_order&(1<<6)) { // r^-6 kspace - if (strcmp(force->kspace_style,"ewald/n") && strcmp(force->kspace_style, "pppm_disp")) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - g_ewald_6 = force->kspace->g_ewald_6; - } + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); + if (ewald_order&(1<<1)) g_ewald = force->kspace->g_ewald; + if (ewald_order&(1<<6)) g_ewald_6 = force->kspace->g_ewald_6; // setup force tables diff --git a/src/KSPACE/pair_coul_long.cpp b/src/KSPACE/pair_coul_long.cpp index a9267cf48a..1aa34041d5 100644 --- a/src/KSPACE/pair_coul_long.cpp +++ b/src/KSPACE/pair_coul_long.cpp @@ -46,6 +46,7 @@ using namespace LAMMPS_NS; PairCoulLong::PairCoulLong(LAMMPS *lmp) : Pair(lmp) { + ewaldflag = pppmflag = 1; ftable = NULL; } @@ -251,7 +252,7 @@ void PairCoulLong::init_style() // insure use of KSpace long-range solver, set g_ewald if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); + error->all(FLERR,"Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables diff --git a/src/KSPACE/pair_coul_msm.cpp b/src/KSPACE/pair_coul_msm.cpp index 46742252f1..7edf074d12 100644 --- a/src/KSPACE/pair_coul_msm.cpp +++ b/src/KSPACE/pair_coul_msm.cpp @@ -38,7 +38,8 @@ using namespace LAMMPS_NS; PairCoulMSM::PairCoulMSM(LAMMPS *lmp) : PairCoulLong(lmp) { - + ewaldflag = pppmflag = 0; + msmflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_charmm_coul_long.cpp b/src/KSPACE/pair_lj_charmm_coul_long.cpp index 6f5b589d71..e81d4784f8 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_long.cpp @@ -48,6 +48,7 @@ using namespace LAMMPS_NS; PairLJCharmmCoulLong::PairLJCharmmCoulLong(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; + ewaldflag = pppmflag = 1; ftable = NULL; implicit = 0; mix_flag = ARITHMETIC; @@ -703,7 +704,8 @@ void PairLJCharmmCoulLong::coeff(int narg, char **arg) void PairLJCharmmCoulLong::init_style() { if (!atom->q_flag) - error->all(FLERR,"Pair style lj/charmm/coul/long requires atom attribute q"); + error->all(FLERR, + "Pair style lj/charmm/coul/long requires atom attribute q"); // request regular or rRESPA neighbor lists @@ -768,7 +770,7 @@ void PairLJCharmmCoulLong::init_style() // insure use of KSpace long-range solver, set g_ewald if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); + error->all(FLERR,"Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.cpp b/src/KSPACE/pair_lj_charmm_coul_msm.cpp index 3539ded08a..7e777605eb 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_msm.cpp @@ -37,12 +37,15 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairLJCharmmCoulMSM::PairLJCharmmCoulMSM(LAMMPS *lmp) : PairLJCharmmCoulLong(lmp) +PairLJCharmmCoulMSM::PairLJCharmmCoulMSM(LAMMPS *lmp) : + PairLJCharmmCoulLong(lmp) { - + ewaldflag = pppmflag = 0; + msmflag = 1; } /* ---------------------------------------------------------------------- */ + void PairLJCharmmCoulMSM::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype,itable; diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index e83e50e239..0a46c4577c 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -50,6 +50,7 @@ using namespace MathConst; PairLJCutCoulLong::PairLJCutCoulLong(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; + ewaldflag = pppmflag = 1; ftable = NULL; qdist = 0.0; } @@ -701,7 +702,7 @@ void PairLJCutCoulLong::init_style() // insure use of KSpace long-range solver, set g_ewald if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); + error->all(FLERR,"Pair style requires a KSpace style"); g_ewald = force->kspace->g_ewald; // setup force tables diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index 9f7ee72de8..1442953341 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -53,6 +53,7 @@ PairLJCutCoulLongTIP4P::PairLJCutCoulLongTIP4P(LAMMPS *lmp) : { single_enable = 0; respa_enable = 0; + tip4pflag = 1; nmax = 0; hneigh = NULL; @@ -455,10 +456,6 @@ void PairLJCutCoulLongTIP4P::init_style() if (!atom->q_flag) error->all(FLERR, "Pair style lj/cut/coul/long/tip4p requires atom attribute q"); - if ( (strcmp(force->kspace_style,"pppm/tip4p") != 0) && - (strcmp(force->kspace_style,"pppm/tip4p/omp") != 0) && - (strcmp(force->kspace_style,"pppm/tip4p/proxy") != 0) ) - error->all(FLERR,"Pair style is incompatible with KSpace style"); if (force->bond == NULL) error->all(FLERR,"Must use a bond style with TIP4P potential"); if (force->angle == NULL) diff --git a/src/KSPACE/pair_lj_cut_coul_msm.cpp b/src/KSPACE/pair_lj_cut_coul_msm.cpp index 087a90d384..2456ced1e2 100644 --- a/src/KSPACE/pair_lj_cut_coul_msm.cpp +++ b/src/KSPACE/pair_lj_cut_coul_msm.cpp @@ -41,7 +41,8 @@ using namespace MathConst; PairLJCutCoulMSM::PairLJCutCoulMSM(LAMMPS *lmp) : PairLJCutCoulLong(lmp) { - + ewaldflag = pppmflag = 0; + msmflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_disp_coul_long.cpp b/src/KSPACE/pair_lj_disp_coul_long.cpp index 333e3e9db1..8718564cdb 100644 --- a/src/KSPACE/pair_lj_disp_coul_long.cpp +++ b/src/KSPACE/pair_lj_disp_coul_long.cpp @@ -48,6 +48,7 @@ using namespace LAMMPS_NS; PairLJDispCoulLong::PairLJDispCoulLong(LAMMPS *lmp) : Pair(lmp) { + dispersionflag = ewaldflag = pppmflag = 1; respa_enable = 1; ftable = NULL; qdist = 0.0; @@ -278,18 +279,8 @@ void PairLJDispCoulLong::init_style() // ensure use of KSpace long-range solver, set g_ewald - if (ewald_order&(1<<1)) { // r^-1 kspace - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - for (i=0; style1[i]&&strcmp(force->kspace_style, style1[i]); ++i); - if (!style1[i]) error->all(FLERR,"Pair style is incompatible with KSpace style"); - } - if (ewald_order&(1<<6)) { // r^-6 kspace - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - for (i=0; style6[i]&&strcmp(force->kspace_style, style6[i]); ++i); - if (!style6[i]) error->all(FLERR,"Pair style is incompatible with KSpace style"); - } + if (force->kspace == NULL) + error->all(FLERR,"Pair style requires a KSpace style"); if (force->kspace) g_ewald = force->kspace->g_ewald; if (force->kspace) g_ewald_6 = force->kspace->g_ewald_6; diff --git a/src/KSPACE/pair_lj_disp_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_disp_coul_long_tip4p.cpp index f5e1473768..c7abb6b81f 100755 --- a/src/KSPACE/pair_lj_disp_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_disp_coul_long_tip4p.cpp @@ -489,9 +489,6 @@ void PairLJDispCoulLongTIP4P::init_style() error->all(FLERR,"Pair style lj/coul/tip4p requires newton pair on"); if (!atom->q_flag) error->all(FLERR,"Pair style lj/coul/tip4p requires atom attribute q"); - if ( (strcmp(force->kspace_style,"pppm_disp/tip4p") != 0) && - (strcmp(force->kspace_style,"pppm_disp/tip4p/proxy") != 0) ) - error->all(FLERR,"Pair style is incompatible with KSpace style"); if (force->bond == NULL) error->all(FLERR,"Must use a bond style with TIP4P potential"); if (force->angle == NULL) diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index fedc1adfde..174b00f3b7 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -60,7 +60,8 @@ using namespace MathConst; PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command"); - + + pppmflag = 1; group_group_enable = 1; accuracy_relative = atof(arg[0]); @@ -194,9 +195,9 @@ void PPPM::init() scale = 1.0; - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - int itmp=0; + pair_check(); + + int itmp = 0; double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp); if (p_cutoff == NULL) error->all(FLERR,"KSpace style is incompatible with Pair style"); @@ -207,10 +208,7 @@ void PPPM::init() qdist = 0.0; - if ((strcmp(force->kspace_style,"pppm/tip4p") == 0) || - (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0)) { - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (tip4pflag) { double *p_qdist = (double *) force->pair->extract("qdist",itmp); int *p_typeO = (int *) force->pair->extract("typeO",itmp); int *p_typeH = (int *) force->pair->extract("typeH",itmp); @@ -237,16 +235,6 @@ void PPPM::init() alpha = qdist / (cos(0.5*theta) * blen); } - // if we have a /proxy pppm version check if the pair style is compatible - - if ((strcmp(force->kspace_style,"pppm/proxy") == 0) || - (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) { - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (strstr(force->pair_style,"pppm/") == NULL ) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - } - // compute qsum & qsqsum and warn if not charge-neutral qsum = qsqsum = 0.0; diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 01e2655623..bd07be353f 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -63,8 +63,8 @@ PPPM_disp::PPPM_disp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm_disp command"); + pppmflag = dispersionflag = 1; accuracy_relative = atof(arg[0]); - nfactors = 3; factors = new int[nfactors]; @@ -81,7 +81,6 @@ PPPM_disp::PPPM_disp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) csumi = NULL; peratom_allocate_flag = 0; - density_brick = vdx_brick = vdy_brick = vdz_brick = NULL; density_fft = NULL; u_brick = v0_brick = v1_brick = v2_brick = v3_brick = v4_brick = v5_brick = NULL; @@ -118,7 +117,6 @@ PPPM_disp::PPPM_disp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) density_fft_a6 = NULL; u_brick_a6 = v0_brick_a6 = v1_brick_a6 = v2_brick_a6 = v3_brick_a6 = v4_brick_a6 = v5_brick_a6 = NULL; - greensfn = NULL; greensfn_6 = NULL; work1 = work2 = NULL; @@ -137,7 +135,6 @@ PPPM_disp::PPPM_disp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) sf_precoeff1 = sf_precoeff2 = sf_precoeff3 = sf_precoeff4 = sf_precoeff5 = sf_precoeff6 = NULL; sf_precoeff1_6 = sf_precoeff2_6 = sf_precoeff3_6 = sf_precoeff4_6 = sf_precoeff5_6 = sf_precoeff6_6 = NULL; - gf_b = NULL; gf_b_6 = NULL; rho1d = rho_coeff = NULL; @@ -161,7 +158,6 @@ PPPM_disp::PPPM_disp(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) com_order = NULL; split_1 = NULL; split_2 = NULL; - } /* ---------------------------------------------------------------------- @@ -184,7 +180,6 @@ PPPM_disp::~PPPM_disp() memory->destroy(dict_rec); memory->destroy(splitbuf1); memory->destroy(splitbuf2); - } /* ---------------------------------------------------------------------- @@ -230,6 +225,8 @@ void PPPM_disp::init() // Check out whether cutoff and pair style are set + pair_check(); + int tmp; Pair *pair = force->pair; int *ptr = pair ? (int *) pair->extract("ewald_order",tmp) : NULL; @@ -242,7 +239,6 @@ void PPPM_disp::init() double tmp2; MPI_Allreduce(&cutoff, &tmp2,1,MPI_DOUBLE,MPI_SUM,world); - // check out which types of potentials will have to be calculated @@ -303,17 +299,12 @@ void PPPM_disp::init() } } - - // if kspace is TIP4P, extract TIP4P params from pair style // bond/angle are not yet init(), so insure equilibrium request is valid qdist = 0.0; - if ( (strcmp(force->kspace_style,"pppm_disp/tip4p") == 0) || - (strcmp(force->kspace_style,"pppm_tip4p/proxy") == 0) ) { - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (tip4pflag) { int itmp; double *p_qdist = (double *) force->pair->extract("qdist",itmp); int *p_typeO = (int *) force->pair->extract("typeO",itmp); @@ -343,15 +334,16 @@ void PPPM_disp::init() // initialize the pair style to get the coefficients + pair->init(); init_coeffs(); //if g_ewald and g_ewald_6 have not been specified, set some initial value // to avoid problems when calculating the energies! + if (!gewaldflag) g_ewald = 1; if (!gewaldflag_6) g_ewald_6 = 1; - // set accuracy (force units) from accuracy_relative or accuracy_absolute if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute; diff --git a/src/KSPACE/pppm_disp_tip4p.cpp b/src/KSPACE/pppm_disp_tip4p.cpp index e9b5016dca..29ac1dcea2 100755 --- a/src/KSPACE/pppm_disp_tip4p.cpp +++ b/src/KSPACE/pppm_disp_tip4p.cpp @@ -42,7 +42,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ PPPMDISPTIP4P::PPPMDISPTIP4P(LAMMPS *lmp, int narg, char **arg) : - PPPM_disp(lmp, narg, arg) {} + PPPM_disp(lmp, narg, arg) +{ + tip4pflag = 1; +} /* ---------------------------------------------------------------------- */ @@ -51,7 +54,7 @@ void PPPMDISPTIP4P::init() // TIP4P PPPM requires newton on, b/c it computes forces on ghost atoms if (force->newton == 0) - error->all(FLERR,"Kspace style pppm_disp/tip4p requires newton on"); + error->all(FLERR,"Kspace style pppm/disp/tip4p requires newton on"); PPPM_disp::init(); } diff --git a/src/KSPACE/pppm_old.cpp b/src/KSPACE/pppm_old.cpp index 796a053885..6f98a4c2ab 100644 --- a/src/KSPACE/pppm_old.cpp +++ b/src/KSPACE/pppm_old.cpp @@ -60,6 +60,7 @@ PPPMOld::PPPMOld(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) { if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command"); + pppmflag = 1; group_group_enable = 1; accuracy_relative = atof(arg[0]); @@ -155,8 +156,8 @@ void PPPMOld::init() scale = 1.0; - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); + pair_check(); + int itmp=0; double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp); if (p_cutoff == NULL) @@ -168,10 +169,7 @@ void PPPMOld::init() qdist = 0.0; - if ( (strcmp(force->kspace_style,"pppm/tip4p") == 0) || - (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) { - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); + if (tip4pflag) { double *p_qdist = (double *) force->pair->extract("qdist",itmp); int *p_typeO = (int *) force->pair->extract("typeO",itmp); int *p_typeH = (int *) force->pair->extract("typeH",itmp); @@ -198,16 +196,6 @@ void PPPMOld::init() alpha = qdist / (cos(0.5*theta) * blen); } - // if we have a /proxy pppm version check if the pair style is compatible - - if ((strcmp(force->kspace_style,"pppm/proxy") == 0) || - (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) { - if (force->pair == NULL) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - if (strstr(force->pair_style,"pppm/") == NULL ) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - } - // compute qsum & qsqsum and warn if not charge-neutral qsum = qsqsum = 0.0; diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp index 1146a64d27..d6c1ddaf0c 100644 --- a/src/KSPACE/pppm_tip4p.cpp +++ b/src/KSPACE/pppm_tip4p.cpp @@ -40,7 +40,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp, int narg, char **arg) : - PPPM(lmp, narg, arg) {} + PPPM(lmp, narg, arg) +{ + tip4pflag = 1; +} /* ---------------------------------------------------------------------- */