jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

This commit is contained in:
sjplimp 2012-10-02 20:11:55 +00:00
parent 27809e8a24
commit 9f8c6721d0
22 changed files with 143 additions and 161 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/KSPACE/ewald.cpp
View File

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

92
src/KSPACE/ewald_disp.cpp
View File

@ -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 (; zz<zn; --zx, ++zy, ++zz) { // set up z[k]=e^(ik.r)
for (; zz<zn; --zx, ++zy, ++zz) { // set up z[k]=e^(ik.r)
C_RMULT(zy->x, 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; k<nk; ++k) { // compute rho(k)
if (ky!=k->y) { // 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 (; f<fn; f+=3) { // -i*qj*fac*
k = kvec; // Sum[conj(d)-d]
for (; f<fn; f+=3) { // -i*qj*fac*
k = kvec; // Sum[conj(d)-d]
kx = ky = -1; // d = k*conj(ekj)*ek
ke = kenergy;
cek = cek_global;
@ -756,7 +752,7 @@ void EwaldDisp::compute_force()
mu++;
}
for (nh = (h = hvec)+nkvec; h<nh; ++h, ++k) {
if (ky!=k->y) { // 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; k<nkvec; ++k) { // sum over k vectors
for (int k=0; k<nkvec; ++k) { // sum over k vectors
if (func[0]) { // 1/r
sum[0] += *(ke++)*(cek->re*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; h<nh; ++h, ++k) {
if (ky!=k->y) { // 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; k<nkvec; ++k) { // sum over k vectors
for (int k=0; k<nkvec; ++k) { // sum over k vectors
if (func[0]) { // 1/r
register double r = cek->re*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; h<nh; ++h, ++k) {
if (ky!=k->y) { // 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;

17
src/KSPACE/msm.cpp
View File

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

7
src/KSPACE/pair_born_coul_long.cpp
View File

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

6
src/KSPACE/pair_born_coul_msm.cpp
View File

@ -35,7 +35,11 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp) {}
PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp)
{
ewaldflag = pppmflag = 0;
msmflag = 1;
}
/* ---------------------------------------------------------------------- */

10
src/KSPACE/pair_buck_coul_long.cpp
View File

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

6
src/KSPACE/pair_buck_coul_msm.cpp
View File

@ -32,7 +32,11 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp) {}
PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp)
{
ewaldflag = pppmflag = 0;
msmflag = 1;
}
/* ---------------------------------------------------------------------- */

37
src/KSPACE/pair_buck_disp_coul_long.cpp
View File

@ -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<<order; break; // set kspace r^-order
case 2: ewald_off |= 1<<order; // turn r^-order off
case 0: ewald_order |= 1<<order; break;
case 2: ewald_off |= 1<<order;
case 1: break;
}
}
/* ---------------------------------------------------------------------- */
void PairBuckDispCoulLong::settings(int narg, char **arg)
{
@ -86,10 +90,13 @@ void PairBuckDispCoulLong::settings(int narg, char **arg)
ewald_order = 0;
ewald_off = 0;
options(arg, 6);
options(++arg, 1);
if (!comm->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

3
src/KSPACE/pair_coul_long.cpp
View File

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

3
src/KSPACE/pair_coul_msm.cpp
View File

@ -38,7 +38,8 @@ using namespace LAMMPS_NS;
PairCoulMSM::PairCoulMSM(LAMMPS *lmp) : PairCoulLong(lmp)
{
ewaldflag = pppmflag = 0;
msmflag = 1;
}
/* ---------------------------------------------------------------------- */

6
src/KSPACE/pair_lj_charmm_coul_long.cpp
View File

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

7
src/KSPACE/pair_lj_charmm_coul_msm.cpp
View File

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

3
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

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

5
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
View File

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

3
src/KSPACE/pair_lj_cut_coul_msm.cpp
View File

@ -41,7 +41,8 @@ using namespace MathConst;
PairLJCutCoulMSM::PairLJCutCoulMSM(LAMMPS *lmp) : PairLJCutCoulLong(lmp)
{
ewaldflag = pppmflag = 0;
msmflag = 1;
}
/* ---------------------------------------------------------------------- */

15
src/KSPACE/pair_lj_disp_coul_long.cpp
View File

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

3
src/KSPACE/pair_lj_disp_coul_long_tip4p.cpp
View File

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

24
src/KSPACE/pppm.cpp
View File

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

20
src/KSPACE/pppm_disp.cpp
View File

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

7
src/KSPACE/pppm_disp_tip4p.cpp
View File

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

20
src/KSPACE/pppm_old.cpp
View File

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

5
src/KSPACE/pppm_tip4p.cpp
View File

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