diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 821c2191bf..9d4dab9299 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -30,6 +30,8 @@ using namespace LAMMPS_NS; +#define SMALL 0.00001 + #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) @@ -78,8 +80,12 @@ void Ewald::init() // error check + if (domain->triclinic) error->all("Cannot use Ewald with triclinic box"); if (force->dimension == 2) error->all("Cannot use Ewald with 2d simulation"); + if (atom->q == NULL) + error->all("Must use charged atom style with kspace style"); + if (slabflag == 0 && domain->nonperiodic > 0) error->all("Cannot use nonperiodic boundaries with Ewald"); if (slabflag == 1) { @@ -99,20 +105,28 @@ void Ewald::init() double cutoff; pair->extract_long(&cutoff); - // compute qsum & qsqsum + // compute qsum & qsqsum and warn if not charge-neutral + + qsum = qsqsum = 0.0; + for (int i = 0; i < atom->nlocal; i++) { + qsum += atom->q[i]; + qsqsum += atom->q[i]*atom->q[i]; + } double tmp; - - qsum = 0.0; - for (int i = 0; i < atom->nlocal; i++) qsum += atom->q[i]; MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsum = tmp; - - qsqsum = 0.0; - for (int i = 0; i < atom->nlocal; i++) qsqsum += atom->q[i]*atom->q[i]; MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsqsum = tmp; + if (qsqsum == 0.0) + error->all("Cannot use kspace solver on system with no charge"); + if (fabs(qsum) > SMALL && comm->me == 0) { + char str[128]; + sprintf(str,"System is not charge neutral, net charge = %g",qsum); + error->warning(str); + } + // setup K-space resolution g_ewald = (1.35 - 0.15*log(precision))/cutoff; diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index 3c37218a92..540f9d5d64 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -250,17 +250,17 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) delx1 = x[i][0] - x2[0]; dely1 = x[i][1] - x2[1]; delz1 = x[i][2] - x2[2]; - domain->minimum_image(&delx1,&dely1,&delz1); + domain->minimum_image(delx1,dely1,delz1); delx2 = x[iH1][0] - x2[0]; dely2 = x[iH1][1] - x2[1]; delz2 = x[iH1][2] - x2[2]; - domain->minimum_image(&delx2,&dely2,&delz2); + domain->minimum_image(delx2,dely2,delz2); delx3 = x[iH2][0] - x2[0]; dely3 = x[iH2][1] - x2[1]; delz3 = x[iH2][2] - x2[2]; - domain->minimum_image(&delx3,&dely3,&delz3); + domain->minimum_image(delx3,dely3,delz3); tvirial[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); tvirial[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); @@ -329,17 +329,17 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) delx1 = x[j][0] - x1[0]; dely1 = x[j][1] - x1[1]; delz1 = x[j][2] - x1[2]; - domain->minimum_image(&delx1,&dely1,&delz1); + domain->minimum_image(delx1,dely1,delz1); delx2 = x[jH1][0] - x1[0]; dely2 = x[jH1][1] - x1[1]; delz2 = x[jH1][2] - x1[2]; - domain->minimum_image(&delx2,&dely2,&delz2); + domain->minimum_image(delx2,dely2,delz2); delx3 = x[jH2][0] - x1[0]; dely3 = x[jH2][1] - x1[1]; delz3 = x[jH2][2] - x1[2]; - domain->minimum_image(&delx3,&dely3,&delz3); + domain->minimum_image(delx3,dely3,delz3); tvirial[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); tvirial[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); @@ -518,12 +518,12 @@ void PairLJCutCoulLongTIP4P::find_M(int i, int &iH1, int &iH2, double *xM) double delx1 = x[iH1][0] - x[i][0]; double dely1 = x[iH1][1] - x[i][1]; double delz1 = x[iH1][2] - x[i][2]; - domain->minimum_image(&delx1,&dely1,&delz1); + domain->minimum_image(delx1,dely1,delz1); double delx2 = x[iH2][0] - x[i][0]; double dely2 = x[iH2][1] - x[i][1]; double delz2 = x[iH2][2] - x[i][2]; - domain->minimum_image(&delx2,&dely2,&delz2); + domain->minimum_image(delx2,dely2,delz2); xM[0] = x[i][0] + alpha * (delx1 + delx2); xM[1] = x[i][1] + alpha * (dely1 + dely2); diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 3bd7857381..79dec55436 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -36,15 +36,15 @@ using namespace LAMMPS_NS; -#define MIN(a,b) ((a) < (b) ? (a) : (b)) -#define MAX(a,b) ((a) > (b) ? (a) : (b)) - #define MAXORDER 7 #define OFFSET 4096 #define SMALL 0.00001 #define LARGE 10000.0 #define EPS_HOC 1.0e-7 +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + /* ---------------------------------------------------------------------- */ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) @@ -105,8 +105,13 @@ void PPPM::init() // error check + if (domain->triclinic) + error->all("Cannot (yet) use PPPM with triclinic box"); if (force->dimension == 2) error->all("Cannot use PPPM with 2d simulation"); + if (atom->q == NULL) + error->all("Must use charged atom style with kspace style"); + if (slabflag == 0 && domain->nonperiodic > 0) error->all("Cannot use nonperiodic boundaries with PPPM"); if (slabflag == 1) { @@ -169,6 +174,8 @@ void PPPM::init() MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world); qsqsum = tmp; + if (qsqsum == 0.0) + error->all("Cannot use kspace solver on system with no charge"); if (fabs(qsum) > SMALL && me == 0) { char str[128]; sprintf(str,"System is not charge neutral, net charge = %g",qsum); @@ -214,36 +221,60 @@ void PPPM::init() // effectively nlo_in,nhi_in + ghost cells // nlo,nhi = global coords of grid pt to "lower left" of smallest/largest // position a particle in my box can be at - // particle position bound = subbox + skin/2.0 + qdist + // dist[3] = particle position bound = subbox + skin/2.0 + qdist // qdist = offset due to TIP4P fictitious charge + // convert to triclinic if necessary // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping // for slab PPPM, assign z grid as if it were not extended - int nlo,nhi; - double cuthalf = 0.5 * neighbor->skin + qdist; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; + triclinic = domain->triclinic; + double *prd,*sublo,*subhi; + + if (triclinic == 0) { + prd = domain->prd; + boxlo = domain->boxlo; + sublo = domain->sublo; + subhi = domain->subhi; + } else { + prd = domain->prd_lamda; + boxlo = domain->boxlo_lamda; + sublo = domain->sublo_lamda; + subhi = domain->subhi_lamda; + } + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; double zprd_slab = zprd*slab_volfactor; - nlo = static_cast ((domain->subxlo-cuthalf-domain->boxxlo) * + double dist[3]; + double cuthalf = 0.5 * neighbor->skin + qdist; + if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf; + else { + dist[0] = cuthalf/domain->prd[0]; + dist[1] = cuthalf/domain->prd[1]; + dist[2] = cuthalf/domain->prd[2]; + } + + int nlo,nhi; + + nlo = static_cast ((sublo[0]-dist[0]-boxlo[0]) * nx_pppm/xprd + shift) - OFFSET; - nhi = static_cast ((domain->subxhi+cuthalf-domain->boxxlo) * + nhi = static_cast ((subhi[0]+dist[0]-boxlo[0]) * nx_pppm/xprd + shift) - OFFSET; nxlo_out = nlo + nlower; nxhi_out = nhi + nupper; - nlo = static_cast ((domain->subylo-cuthalf-domain->boxylo) * + nlo = static_cast ((sublo[1]-dist[1]-boxlo[1]) * ny_pppm/yprd + shift) - OFFSET; - nhi = static_cast ((domain->subyhi+cuthalf-domain->boxylo) * + nhi = static_cast ((subhi[1]+dist[1]-boxlo[1]) * ny_pppm/yprd + shift) - OFFSET; nylo_out = nlo + nlower; nyhi_out = nhi + nupper; - nlo = static_cast ((domain->subzlo-cuthalf-domain->boxzlo) * + nlo = static_cast ((sublo[2]-dist[2]-boxlo[2]) * nz_pppm/zprd_slab + shift) - OFFSET; - nhi = static_cast ((domain->subzhi+cuthalf-domain->boxzlo) * + nhi = static_cast ((subhi[2]+dist[2]-boxlo[2]) * nz_pppm/zprd_slab + shift) - OFFSET; nzlo_out = nlo + nlower; nzhi_out = nhi + nupper; @@ -426,18 +457,19 @@ void PPPM::init() void PPPM::setup() { int i,j,k,l,m,n; + double *prd; // volume-dependent factors + // adjust z dimension for 2d slab PPPM + // z dimension for 3d PPPM is zprd since slab_volfactor = 1.0 - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - - // adjustment of z dimension for 2d slab PPPM - // 3d PPPM just uses zprd since slab_volfactor = 1.0 + if (triclinic == 0) prd = domain->prd; + else prd = domain->prd_lamda; + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; double zprd_slab = zprd*slab_volfactor; - volume = xprd * yprd * zprd_slab; delxinv = nx_pppm/xprd; @@ -579,6 +611,14 @@ void PPPM::compute(int eflag, int vflag) { int i; + // convert atoms from box to lamda coords + + if (triclinic == 0) boxlo = domain->boxlo; + else { + boxlo = domain->boxlo_lamda; + domain->x2lamda(atom->nlocal); + } + // extend size of nlocal-dependent arrays if necessary if (atom->nlocal > nmax) { @@ -641,6 +681,10 @@ void PPPM::compute(int eflag, int vflag) // 2d slab correction if (slabflag) slabcorr(eflag); + + // convert atoms back from lamda to box coords + + if (triclinic) domain->lamda2x(atom->nlocal); } /* ---------------------------------------------------------------------- @@ -781,7 +825,8 @@ void PPPM::set_grid() double q2 = qsqsum / force->dielectric; double natoms = atom->natoms; - // adjustment of z dimension for 2d slab PPPM + // use xprd,yprd,zprd even if triclinic so grid size is the same + // adjust z dimension for 2d slab PPPM // 3d PPPM just uses zprd since slab_volfactor = 1.0 double xprd = domain->xprd; @@ -1387,9 +1432,6 @@ void PPPM::particle_map() double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; int flag = 0; for (int i = 0; i < nlocal; i++) { @@ -1398,9 +1440,9 @@ void PPPM::particle_map() // current particle coord can be outside global and local box // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 - nx = static_cast ((x[i][0]-boxxlo)*delxinv+shift) - OFFSET; - ny = static_cast ((x[i][1]-boxylo)*delyinv+shift) - OFFSET; - nz = static_cast ((x[i][2]-boxzlo)*delzinv+shift) - OFFSET; + nx = static_cast ((x[i][0]-boxlo[0])*delxinv+shift) - OFFSET; + ny = static_cast ((x[i][1]-boxlo[1])*delyinv+shift) - OFFSET; + nz = static_cast ((x[i][2]-boxlo[2])*delzinv+shift) - OFFSET; part2grid[i][0] = nx; part2grid[i][1] = ny; @@ -1443,18 +1485,15 @@ void PPPM::make_rho() double *q = atom->q; double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; for (int i = 0; i < nlocal; i++) { nx = part2grid[i][0]; ny = part2grid[i][1]; nz = part2grid[i][2]; - dx = nx+shiftone - (x[i][0]-boxxlo)*delxinv; - dy = ny+shiftone - (x[i][1]-boxylo)*delyinv; - dz = nz+shiftone - (x[i][2]-boxzlo)*delzinv; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; compute_rho1d(dx,dy,dz); @@ -1614,18 +1653,15 @@ void PPPM::fieldforce() double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; for (i = 0; i < nlocal; i++) { nx = part2grid[i][0]; ny = part2grid[i][1]; nz = part2grid[i][2]; - dx = nx+shiftone - (x[i][0]-boxxlo)*delxinv; - dy = ny+shiftone - (x[i][1]-boxylo)*delyinv; - dz = nz+shiftone - (x[i][2]-boxzlo)*delzinv; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; compute_rho1d(dx,dy,dz); diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index 2acaa9f30b..021d6dfc4f 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -66,6 +66,8 @@ class PPPM : public KSpace { int **part2grid; // storage for particle -> grid mapping int nmax; + int triclinic; // domain settings, orthog or triclinic + double *boxlo; // TIP4P settings int typeH,typeO; // atom types of TIP4P water H and O atoms double qdist; // distance from O site to negative charge diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp index 3486e95287..dfab007945 100644 --- a/src/KSPACE/pppm_tip4p.cpp +++ b/src/KSPACE/pppm_tip4p.cpp @@ -45,9 +45,6 @@ void PPPMTIP4P::particle_map() int *type = atom->type; double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; int flag = 0; for (int i = 0; i < nlocal; i++) { @@ -60,9 +57,9 @@ void PPPMTIP4P::particle_map() // current particle coord can be outside global and local box // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1 - nx = static_cast ((xi[0]-boxxlo)*delxinv+shift) - OFFSET; - ny = static_cast ((xi[1]-boxylo)*delyinv+shift) - OFFSET; - nz = static_cast ((xi[2]-boxzlo)*delzinv+shift) - OFFSET; + nx = static_cast ((xi[0]-boxlo[0])*delxinv+shift) - OFFSET; + ny = static_cast ((xi[1]-boxlo[1])*delyinv+shift) - OFFSET; + nz = static_cast ((xi[2]-boxlo[2])*delzinv+shift) - OFFSET; part2grid[i][0] = nx; part2grid[i][1] = ny; @@ -107,9 +104,6 @@ void PPPMTIP4P::make_rho() double *q = atom->q; double **x = atom->x; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; for (int i = 0; i < nlocal; i++) { if (type[i] == typeO) { @@ -120,9 +114,9 @@ void PPPMTIP4P::make_rho() nx = part2grid[i][0]; ny = part2grid[i][1]; nz = part2grid[i][2]; - dx = nx+shiftone - (xi[0]-boxxlo)*delxinv; - dy = ny+shiftone - (xi[1]-boxylo)*delyinv; - dz = nz+shiftone - (xi[2]-boxzlo)*delzinv; + dx = nx+shiftone - (xi[0]-boxlo[0])*delxinv; + dy = ny+shiftone - (xi[1]-boxlo[1])*delyinv; + dz = nz+shiftone - (xi[2]-boxlo[2])*delzinv; compute_rho1d(dx,dy,dz); @@ -167,9 +161,6 @@ void PPPMTIP4P::fieldforce() double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - double boxxlo = domain->boxxlo; - double boxylo = domain->boxylo; - double boxzlo = domain->boxzlo; for (i = 0; i < nlocal; i++) { if (type[i] == typeO) { @@ -180,9 +171,9 @@ void PPPMTIP4P::fieldforce() nx = part2grid[i][0]; ny = part2grid[i][1]; nz = part2grid[i][2]; - dx = nx+shiftone - (xi[0]-boxxlo)*delxinv; - dy = ny+shiftone - (xi[1]-boxylo)*delyinv; - dz = nz+shiftone - (xi[2]-boxzlo)*delzinv; + dx = nx+shiftone - (xi[0]-boxlo[0])*delxinv; + dy = ny+shiftone - (xi[1]-boxlo[1])*delyinv; + dz = nz+shiftone - (xi[2]-boxlo[2])*delzinv; compute_rho1d(dx,dy,dz); @@ -251,12 +242,12 @@ void PPPMTIP4P::find_M(int i, int &iH1, int &iH2, double *xM) double delx1 = x[iH1][0] - x[i][0]; double dely1 = x[iH1][1] - x[i][1]; double delz1 = x[iH1][2] - x[i][2]; - domain->minimum_image(&delx1,&dely1,&delz1); + domain->minimum_image(delx1,dely1,delz1); double delx2 = x[iH2][0] - x[i][0]; double dely2 = x[iH2][1] - x[i][1]; double delz2 = x[iH2][2] - x[i][2]; - domain->minimum_image(&delx2,&dely2,&delz2); + domain->minimum_image(delx2,dely2,delz2); xM[0] = x[i][0] + alpha * (delx1 + delx2); xM[1] = x[i][1] + alpha * (dely1 + dely2);