diff --git a/src/GPU/pair_lj_cut_coul_msm_gpu.cpp b/src/GPU/pair_lj_cut_coul_msm_gpu.cpp index 78e0243c36..4d9486d1fa 100644 --- a/src/GPU/pair_lj_cut_coul_msm_gpu.cpp +++ b/src/GPU/pair_lj_cut_coul_msm_gpu.cpp @@ -136,6 +136,9 @@ void PairLJCutCoulMSMGPU::init_style() if (force->newton_pair) error->all(FLERR,"Cannot use newton pair with lj/cut/coul/msm/gpu pair style"); + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' with GPU MSM Pair styles"); + // Repeat cutsq calculation because done after call to init_style double maxcut = -1.0; double cut; diff --git a/src/GPU/pair_lj_cut_coul_msm_gpu.h b/src/GPU/pair_lj_cut_coul_msm_gpu.h index f77a5abca7..c125edcd6d 100644 --- a/src/GPU/pair_lj_cut_coul_msm_gpu.h +++ b/src/GPU/pair_lj_cut_coul_msm_gpu.h @@ -56,4 +56,8 @@ E: Cannot use newton pair with lj/cut/coul/msm/gpu pair style Self-explanatory. +E: Must use 'kspace_modify pressure/scalar no' with GPU MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with GPU MSM Pair styles. + */ diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index a99d5bb1fb..809de53d7a 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -100,6 +100,7 @@ MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) levels = 0; peratom_allocate_flag = 0; + scalar_pressure_flag = 1; order = 10; } @@ -416,7 +417,7 @@ void MSM::setup() } else { get_g_direct(); if (domain->nonperiodic) get_g_direct_top(levels-1); - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { get_virial_direct(); if (domain->nonperiodic) get_virial_direct_top(levels-1); } @@ -463,6 +464,19 @@ void MSM::compute(int eflag, int vflag) else evflag = evflag_atom = eflag_global = vflag_global = eflag_atom = vflag_atom = eflag_either = vflag_either = 0; + if (scalar_pressure_flag && vflag_either) { + if (vflag_atom) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' to obtain " + "per-atom virial with kspace_style MSM"); + + // must switch on global energy computation if not already on + + if (eflag == 0 || eflag == 2) { + eflag++; + ev_setup(eflag,vflag); + } + } + // invoke allocate_peratom() if needed for first time if (vflag_atom && !peratom_allocate_flag) { @@ -588,12 +602,17 @@ void MSM::compute(int eflag, int vflag) // total long-range virial - if (vflag_global) { + if (vflag_global && !scalar_pressure_flag) { double virial_all[6]; MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*virial_all[i]; } + // fast compute of scalar pressure (if requested) + + if (scalar_pressure_flag && vflag_global) + for (i = 0; i < 3; i++) virial[i] = energy/3.0; + // per-atom energy/virial // energy includes self-energy correction @@ -1588,7 +1607,7 @@ void MSM::direct(int n) qtmp = qgridn[icz][icy][icx]; // charge on center grid point esum = 0.0; - if (vflag_either) + if (vflag_either && !scalar_pressure_flag) v0sum = v1sum = v2sum = v3sum = v4sum = v5sum = 0.0; // use hemisphere to avoid double computation of pair-wise @@ -1612,7 +1631,7 @@ void MSM::direct(int n) esum += gtmp * qtmp2; ekj[ii] += gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_directn[k] * qtmp2; v1sum += v1_directn[k] * qtmp2; v2sum += v2_directn[k] * qtmp2; @@ -1644,7 +1663,7 @@ void MSM::direct(int n) esum += gtmp * qtmp2; ekj[ii] += gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_directn[k] * qtmp2; v1sum += v1_directn[k] * qtmp2; v2sum += v2_directn[k] * qtmp2; @@ -1675,7 +1694,7 @@ void MSM::direct(int n) esum += gtmp * qtmp2; ekj[ii] += gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_directn[k] * qtmp2; v1sum += v1_directn[k] * qtmp2; v2sum += v2_directn[k] * qtmp2; @@ -1717,7 +1736,7 @@ void MSM::direct(int n) if (evflag) { qtmp = qgridn[icz][icy][icx]; if (eflag_global) energy += 2.0 * esum * qtmp; - if (vflag_global) { + if (vflag_global && !scalar_pressure_flag) { virial[0] += 2.0 * v0sum * qtmp; virial[1] += 2.0 * v1sum * qtmp; virial[2] += 2.0 * v2sum * qtmp; @@ -1947,7 +1966,7 @@ void MSM::direct_top(int n) qtmp = qgridn[icz][icy][icx]; esum = 0.0; - if (vflag_either) + if (vflag_either && !scalar_pressure_flag) v0sum = v1sum = v2sum = v3sum = v4sum = v5sum = 0.0; // use hemisphere to avoid double computation of pair-wise @@ -1971,7 +1990,7 @@ void MSM::direct_top(int n) esum += gtmp * qtmp2; ekj[ii] += gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_direct_top[k] * qtmp2; v1sum += v1_direct_top[k] * qtmp2; v2sum += v2_direct_top[k] * qtmp2; @@ -2003,7 +2022,7 @@ void MSM::direct_top(int n) esum += gtmp * qtmp2; ekj[ii] += gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_direct_top[k] * qtmp2; v1sum += v1_direct_top[k] * qtmp2; v2sum += v2_direct_top[k] * qtmp2; @@ -2034,7 +2053,7 @@ void MSM::direct_top(int n) esum += gtmp * qtmp2; ekj[ii] += gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_direct_top[k] * qtmp2; v1sum += v1_direct_top[k] * qtmp2; v2sum += v2_direct_top[k] * qtmp2; @@ -2057,7 +2076,7 @@ void MSM::direct_top(int n) esum += 0.5 * gtmp * qtmp; egridn[icz][icy][icx] += 0.5 * gtmp * qtmp; - if (vflag_either) { + if (vflag_either && !scalar_pressure_flag) { v0sum += v0_direct_top[k] * qtmp; v1sum += v1_direct_top[k] * qtmp; v2sum += v2_direct_top[k] * qtmp; @@ -2084,7 +2103,7 @@ void MSM::direct_top(int n) if (evflag) { qtmp = qgridn[icz][icy][icx]; if (eflag_global) energy += 2.0 * esum * qtmp; - if (vflag_global) { + if (vflag_global && !scalar_pressure_flag) { virial[0] += 2.0 * v0sum * qtmp; virial[1] += 2.0 * v1sum * qtmp; virial[2] += 2.0 * v2sum * qtmp; diff --git a/src/KSPACE/msm.h b/src/KSPACE/msm.h index 7bee1ebc48..3d333e9073 100644 --- a/src/KSPACE/msm.h +++ b/src/KSPACE/msm.h @@ -217,6 +217,11 @@ The global MSM grid is larger than OFFSET in one or more dimensions. OFFSET is currently set to 16384. You likely need to decrease the requested accuracy. +E: Must use 'kspace_modify pressure/scalar no' to obtain per-atom virial +with kspace_style MSM + +The kspace scalar pressure option cannot be used to obtain per-atom virial. + E: Out of range atoms - cannot compute MSM One or more atoms are attempting to map their charge to a MSM grid point diff --git a/src/KSPACE/msm_cg.cpp b/src/KSPACE/msm_cg.cpp index 527000f651..23d719ef04 100644 --- a/src/KSPACE/msm_cg.cpp +++ b/src/KSPACE/msm_cg.cpp @@ -72,6 +72,10 @@ MSMCG::~MSMCG() void MSMCG::compute(int eflag, int vflag) { + if (scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' with " + "kspace_style msm/cg"); + const double * const q = atom->q; const int nlocal = atom->nlocal; int i,j,n; diff --git a/src/KSPACE/msm_cg.h b/src/KSPACE/msm_cg.h index a1009ce513..44c5fd3704 100644 --- a/src/KSPACE/msm_cg.h +++ b/src/KSPACE/msm_cg.h @@ -56,6 +56,10 @@ Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. +E: Must use 'kspace_modify pressure/scalar no' with kspace_style msm/cg + +The kspace scalar pressure option is not compatible with kspace_style msm/cg. + E: Out of range atoms - cannot compute MSM One or more atoms are attempting to map their charge to a MSM grid point diff --git a/src/KSPACE/pair_born_coul_msm.cpp b/src/KSPACE/pair_born_coul_msm.cpp index a94b753bfd..cbb11cdcfd 100644 --- a/src/KSPACE/pair_born_coul_msm.cpp +++ b/src/KSPACE/pair_born_coul_msm.cpp @@ -39,6 +39,15 @@ PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp) { ewaldflag = pppmflag = 0; msmflag = 1; + nmax = 0; + ftmp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairBornCoulMSM::~PairBornCoulMSM() +{ + if (ftmp) memory->destroy(ftmp); } /* ---------------------------------------------------------------------- */ @@ -46,11 +55,31 @@ PairBornCoulMSM::PairBornCoulMSM(LAMMPS *lmp) : PairBornCoulLong(lmp) void PairBornCoulMSM::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair,fcoul; double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj; double egamma,fgamma,prefactor; double r,rexp; int *ilist,*jlist,*numneigh,**firstneigh; + int eflag_old = eflag; + + if (force->kspace->scalar_pressure_flag && vflag) { + if (vflag > 2) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' to " + "obtain per-atom virial with kspace_style MSM"); + + if (atom->nmax > nmax) { + if (ftmp) memory->destroy(ftmp); + nmax = atom->nmax; + memory->create(ftmp,nmax,3,"pair:ftmp"); + } + memset(&ftmp[0][0],0,nmax*3*sizeof(double)); + + // must switch on global energy computation if not already on + + if (eflag == 0 || eflag == 2) { + eflag++; + } + } evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -114,15 +143,41 @@ void PairBornCoulMSM::compute(int eflag, int vflag) + born3[itype][jtype]*r2inv*r6inv; } else forceborn = 0.0; - fpair = (forcecoul + factor_lj*forceborn) * r2inv; + if (!(force->kspace->scalar_pressure_flag && vflag)) { + fpair = (forcecoul + factor_lj*forceborn) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + } else { + // separate Born and Coulombic forces + + fpair = (factor_lj*forceborn) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + fcoul = (forcecoul) * r2inv; + + ftmp[i][0] += delx*fcoul; + ftmp[i][1] += dely*fcoul; + ftmp[i][2] += delz*fcoul; + if (newton_pair || j < nlocal) { + ftmp[j][0] -= delx*fcoul; + ftmp[j][1] -= dely*fcoul; + ftmp[j][2] -= delz*fcoul; + } } if (eflag) { @@ -130,7 +185,7 @@ void PairBornCoulMSM::compute(int eflag, int vflag) ecoul = prefactor*egamma; if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; - if (rsq < cut_ljsq[itype][jtype]) { + if (eflag_old && rsq < cut_ljsq[itype][jtype]) { evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv + d[itype][jtype]*r6inv*r2inv - offset[itype][jtype]; evdwl *= factor_lj; @@ -144,6 +199,15 @@ void PairBornCoulMSM::compute(int eflag, int vflag) } if (vflag_fdotr) virial_fdotr_compute(); + + if (force->kspace->scalar_pressure_flag && vflag) { + for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0; + for (int i = 0; i < nmax; i++) { + f[i][0] += ftmp[i][0]; + f[i][1] += ftmp[i][1]; + f[i][2] += ftmp[i][2]; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_born_coul_msm.h b/src/KSPACE/pair_born_coul_msm.h index b4d8368d1f..2828b49264 100644 --- a/src/KSPACE/pair_born_coul_msm.h +++ b/src/KSPACE/pair_born_coul_msm.h @@ -27,10 +27,13 @@ namespace LAMMPS_NS { class PairBornCoulMSM : public PairBornCoulLong { public: PairBornCoulMSM(class LAMMPS *); - virtual ~PairBornCoulMSM(){}; + virtual ~PairBornCoulMSM(); virtual void compute(int, int); virtual double single(int, int, int, int, double, double, double, double &); virtual void *extract(const char *, int &); + protected: + int nmax; + double **ftmp; }; @@ -65,4 +68,8 @@ E: Pair style is incompatible with KSpace style If a pair style with a long-range Coulombic component is selected, then a kspace style must also be used. +E: Must use 'kspace_modify pressure/scalar no' to obtain per-atom virial with kspace_style MSM + +The kspace scalar pressure option cannot be used to obtain per-atom virial. + */ diff --git a/src/KSPACE/pair_buck_coul_msm.cpp b/src/KSPACE/pair_buck_coul_msm.cpp index 4a8febae5a..efd0bc8b97 100644 --- a/src/KSPACE/pair_buck_coul_msm.cpp +++ b/src/KSPACE/pair_buck_coul_msm.cpp @@ -36,6 +36,15 @@ PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp) { ewaldflag = pppmflag = 0; msmflag = 1; + nmax = 0; + ftmp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairBuckCoulMSM::~PairBuckCoulMSM() +{ + if (ftmp) memory->destroy(ftmp); } /* ---------------------------------------------------------------------- */ @@ -43,11 +52,31 @@ PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp) void PairBuckCoulMSM::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair,fcoul; double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj; double egamma,fgamma,prefactor; double r,rexp; int *ilist,*jlist,*numneigh,**firstneigh; + int eflag_old = eflag; + + if (force->kspace->scalar_pressure_flag && vflag) { + if (vflag > 2) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "to obtain per-atom virial with kspace_style MSM"); + + if (atom->nmax > nmax) { + if (ftmp) memory->destroy(ftmp); + nmax = atom->nmax; + memory->create(ftmp,nmax,3,"pair:ftmp"); + } + memset(&ftmp[0][0],0,nmax*3*sizeof(double)); + + // must switch on global energy computation if not already on + + if (eflag == 0 || eflag == 2) { + eflag++; + } + } evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -110,15 +139,41 @@ void PairBuckCoulMSM::compute(int eflag, int vflag) forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv; } else forcebuck = 0.0; - fpair = (forcecoul + factor_lj*forcebuck) * r2inv; + if (!(force->kspace->scalar_pressure_flag && vflag)) { + fpair = (forcecoul + factor_lj*forcebuck) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + } else { + // separate Buck and Coulombic forces + + fpair = (factor_lj*forcebuck) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + fcoul = (forcecoul) * r2inv; + + ftmp[i][0] += delx*fcoul; + ftmp[i][1] += dely*fcoul; + ftmp[i][2] += delz*fcoul; + if (newton_pair || j < nlocal) { + ftmp[j][0] -= delx*fcoul; + ftmp[j][1] -= dely*fcoul; + ftmp[j][2] -= delz*fcoul; + } } if (eflag) { @@ -126,7 +181,7 @@ void PairBuckCoulMSM::compute(int eflag, int vflag) ecoul = prefactor*egamma; if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; - if (rsq < cut_ljsq[itype][jtype]) { + if (eflag_old && rsq < cut_ljsq[itype][jtype]) { evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv - offset[itype][jtype]; evdwl *= factor_lj; @@ -140,6 +195,15 @@ void PairBuckCoulMSM::compute(int eflag, int vflag) } if (vflag_fdotr) virial_fdotr_compute(); + + if (force->kspace->scalar_pressure_flag && vflag) { + for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0; + for (int i = 0; i < nmax; i++) { + f[i][0] += ftmp[i][0]; + f[i][1] += ftmp[i][1]; + f[i][2] += ftmp[i][2]; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_buck_coul_msm.h b/src/KSPACE/pair_buck_coul_msm.h index 95f63268a8..5bce892172 100644 --- a/src/KSPACE/pair_buck_coul_msm.h +++ b/src/KSPACE/pair_buck_coul_msm.h @@ -27,11 +27,13 @@ namespace LAMMPS_NS { class PairBuckCoulMSM : public PairBuckCoulLong { public: PairBuckCoulMSM(class LAMMPS *); - virtual ~PairBuckCoulMSM(){}; + virtual ~PairBuckCoulMSM(); virtual void compute(int, int); virtual double single(int, int, int, int, double, double, double, double &); virtual void *extract(const char *, int &); - + protected: + int nmax; + double **ftmp; }; } @@ -65,4 +67,8 @@ E: Pair style is incompatible with KSpace style If a pair style with a long-range Coulombic component is selected, then a kspace style must also be used. +E: Must use 'kspace_modify pressure/scalar no' to obtain per-atom virial with kspace_style MSM + +The kspace scalar pressure option cannot be used to obtain per-atom virial. + */ diff --git a/src/KSPACE/pair_coul_msm.cpp b/src/KSPACE/pair_coul_msm.cpp index b390fb3e15..e27714f59d 100644 --- a/src/KSPACE/pair_coul_msm.cpp +++ b/src/KSPACE/pair_coul_msm.cpp @@ -53,7 +53,19 @@ void PairCoulMSM::compute(int eflag, int vflag) double egamma,fgamma,prefactor; int *ilist,*jlist,*numneigh,**firstneigh; double rsq; - + + if (force->kspace->scalar_pressure_flag && vflag) { + if (vflag > 2) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "to obtain per-atom virial with kspace_style MSM"); + + // must switch on global energy computation if not already on + + if (eflag == 0 || eflag == 2) { + eflag++; + } + } + ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; @@ -140,13 +152,19 @@ void PairCoulMSM::compute(int eflag, int vflag) if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } + if (force->kspace->scalar_pressure_flag) + fpair = 0.0; + if (evflag) ev_tally(i,j,nlocal,newton_pair, 0.0,ecoul,fpair,delx,dely,delz); } } } - if (vflag_fdotr) virial_fdotr_compute(); + if (vflag_fdotr && !force->kspace->scalar_pressure_flag) + virial_fdotr_compute(); + if (force->kspace->scalar_pressure_flag && vflag) + for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0; } /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_coul_msm.h b/src/KSPACE/pair_coul_msm.h index 1dad2b4852..aaa3a63d5d 100644 --- a/src/KSPACE/pair_coul_msm.h +++ b/src/KSPACE/pair_coul_msm.h @@ -64,4 +64,8 @@ E: Pair style is incompatible with KSpace style If a pair style with a long-range Coulombic component is selected, then a kspace style must also be used. +E: Must use 'kspace_modify pressure/scalar no' to obtain per-atom virial with kspace_style MSM + +The kspace scalar pressure option cannot be used to obtain per-atom virial. + */ diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.cpp b/src/KSPACE/pair_lj_charmm_coul_msm.cpp index 1d3084764e..c8f0841608 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_msm.cpp @@ -42,6 +42,15 @@ PairLJCharmmCoulMSM::PairLJCharmmCoulMSM(LAMMPS *lmp) : { ewaldflag = pppmflag = 0; msmflag = 1; + nmax = 0; + ftmp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCharmmCoulMSM::~PairLJCharmmCoulMSM() +{ + if (ftmp) memory->destroy(ftmp); } /* ---------------------------------------------------------------------- */ @@ -49,13 +58,33 @@ PairLJCharmmCoulMSM::PairLJCharmmCoulMSM(LAMMPS *lmp) : void PairLJCharmmCoulMSM::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair,fcoul; double fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; double egamma,fgamma,prefactor; double philj,switch1,switch2; int *ilist,*jlist,*numneigh,**firstneigh; double rsq; + int eflag_old = eflag; + + if (force->kspace->scalar_pressure_flag && vflag) { + if (vflag > 2) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "to obtain per-atom virial with kspace_style MSM"); + + if (atom->nmax > nmax) { + if (ftmp) memory->destroy(ftmp); + nmax = atom->nmax; + memory->create(ftmp,nmax,3,"pair:ftmp"); + } + memset(&ftmp[0][0],0,nmax*3*sizeof(double)); + + // must switch on global energy computation if not already on + + if (eflag == 0 || eflag == 2) { + eflag++; + } + } evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -140,15 +169,41 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag) } } else forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + if (!(force->kspace->scalar_pressure_flag && vflag)) { + fpair = (forcecoul + factor_lj*forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + } else { + // separate LJ and Coulombic forces + + fpair = (factor_lj*forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + fcoul = (forcecoul) * r2inv; + + ftmp[i][0] += delx*fcoul; + ftmp[i][1] += dely*fcoul; + ftmp[i][2] += delz*fcoul; + if (newton_pair || j < nlocal) { + ftmp[j][0] -= delx*fcoul; + ftmp[j][1] -= dely*fcoul; + ftmp[j][2] -= delz*fcoul; + } } if (eflag) { @@ -162,7 +217,7 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag) if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; - if (rsq < cut_ljsq) { + if (eflag_old && rsq < cut_ljsq) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); if (rsq > cut_lj_innersq) { switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) * @@ -180,6 +235,15 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag) } if (vflag_fdotr) virial_fdotr_compute(); + + if (force->kspace->scalar_pressure_flag && vflag) { + for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0; + for (int i = 0; i < nmax; i++) { + f[i][0] += ftmp[i][0]; + f[i][1] += ftmp[i][1]; + f[i][2] += ftmp[i][2]; + } + } } /* ---------------------------------------------------------------------- */ @@ -196,6 +260,10 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) int *ilist,*jlist,*numneigh,**firstneigh; double rsq; + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "for rRESPA with kspace_style MSM"); + evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.h b/src/KSPACE/pair_lj_charmm_coul_msm.h index 2efe34777b..e9dac633f2 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.h +++ b/src/KSPACE/pair_lj_charmm_coul_msm.h @@ -27,11 +27,14 @@ namespace LAMMPS_NS { class PairLJCharmmCoulMSM : public PairLJCharmmCoulLong { public: PairLJCharmmCoulMSM(class LAMMPS *); - virtual ~PairLJCharmmCoulMSM(){}; + virtual ~PairLJCharmmCoulMSM(); virtual void compute(int, int); virtual double single(int, int, int, int, double, double, double, double &); virtual void compute_outer(int, int); virtual void *extract(const char *, int &); + protected: + int nmax; + double **ftmp; }; } @@ -59,6 +62,15 @@ E: Pair inner cutoff >= Pair outer cutoff The specified cutoffs for the pair style are inconsistent. +E: Must use 'kspace_modify pressure/scalar no' to obtain per-atom virial +with kspace_style MSM + +The kspace scalar pressure option cannot be used to obtain per-atom virial. + +E: Must use 'kspace_modify pressure/scalar no' for rRESPA with kspace_style MSM + +The kspace scalar pressure option cannot (yet) be used with rRESPA. + E: Pair cutoff < Respa interior cutoff One or more pairwise cutoffs are too short to use with the specified diff --git a/src/KSPACE/pair_lj_cut_coul_msm.cpp b/src/KSPACE/pair_lj_cut_coul_msm.cpp index 84d764c14b..d075b7db6d 100644 --- a/src/KSPACE/pair_lj_cut_coul_msm.cpp +++ b/src/KSPACE/pair_lj_cut_coul_msm.cpp @@ -43,6 +43,15 @@ PairLJCutCoulMSM::PairLJCutCoulMSM(LAMMPS *lmp) : PairLJCutCoulLong(lmp) { ewaldflag = pppmflag = 0; msmflag = 1; + nmax = 0; + ftmp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairLJCutCoulMSM::~PairLJCutCoulMSM() +{ + if (ftmp) memory->destroy(ftmp); } /* ---------------------------------------------------------------------- */ @@ -50,12 +59,32 @@ PairLJCutCoulMSM::PairLJCutCoulMSM(LAMMPS *lmp) : PairLJCutCoulLong(lmp) void PairLJCutCoulMSM::compute(int eflag, int vflag) { int i,ii,j,jj,inum,jnum,itype,jtype,itable; - double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair,fcoul; double fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; double egamma,fgamma,prefactor; int *ilist,*jlist,*numneigh,**firstneigh; double rsq; + int eflag_old = eflag; + + if (force->kspace->scalar_pressure_flag && vflag) { + if (vflag > 2) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "to obtain per-atom virial with kspace_style MSM"); + + if (atom->nmax > nmax) { + if (ftmp) memory->destroy(ftmp); + nmax = atom->nmax; + memory->create(ftmp,nmax,3,"pair:ftmp"); + } + memset(&ftmp[0][0],0,nmax*3*sizeof(double)); + + // must switch on global energy computation if not already on + + if (eflag == 0 || eflag == 2) { + eflag++; + } + } evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); @@ -132,15 +161,41 @@ void PairLJCutCoulMSM::compute(int eflag, int vflag) forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); } else forcelj = 0.0; - fpair = (forcecoul + factor_lj*forcelj) * r2inv; + if (!(force->kspace->scalar_pressure_flag && vflag)) { + fpair = (forcecoul + factor_lj*forcelj) * r2inv; - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + } else { + // separate LJ and Coulombic forces + + fpair = (factor_lj*forcelj) * r2inv; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + fcoul = (forcecoul) * r2inv; + + ftmp[i][0] += delx*fcoul; + ftmp[i][1] += dely*fcoul; + ftmp[i][2] += delz*fcoul; + if (newton_pair || j < nlocal) { + ftmp[j][0] -= delx*fcoul; + ftmp[j][1] -= dely*fcoul; + ftmp[j][2] -= delz*fcoul; + } } if (eflag) { @@ -154,7 +209,7 @@ void PairLJCutCoulMSM::compute(int eflag, int vflag) if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; - if (rsq < cut_ljsq[itype][jtype]) { + if (eflag_old && rsq < cut_ljsq[itype][jtype]) { evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; @@ -168,6 +223,15 @@ void PairLJCutCoulMSM::compute(int eflag, int vflag) } if (vflag_fdotr) virial_fdotr_compute(); + + if (force->kspace->scalar_pressure_flag && vflag) { + for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0; + for (int i = 0; i < nmax; i++) { + f[i][0] += ftmp[i][0]; + f[i][1] += ftmp[i][1]; + f[i][2] += ftmp[i][2]; + } + } } /* ---------------------------------------------------------------------- */ @@ -183,6 +247,10 @@ void PairLJCutCoulMSM::compute_outer(int eflag, int vflag) int *ilist,*jlist,*numneigh,**firstneigh; double rsq; + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "for rRESPA with kspace_style MSM"); + evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; diff --git a/src/KSPACE/pair_lj_cut_coul_msm.h b/src/KSPACE/pair_lj_cut_coul_msm.h index 3c88e22896..b1e9fe1c0b 100644 --- a/src/KSPACE/pair_lj_cut_coul_msm.h +++ b/src/KSPACE/pair_lj_cut_coul_msm.h @@ -27,11 +27,14 @@ namespace LAMMPS_NS { class PairLJCutCoulMSM : public PairLJCutCoulLong { public: PairLJCutCoulMSM(class LAMMPS *); - virtual ~PairLJCutCoulMSM(){}; + virtual ~PairLJCutCoulMSM(); virtual void compute(int, int); virtual double single(int, int, int, int, double, double, double, double &); virtual void compute_outer(int, int); virtual void *extract(const char *, int &); + protected: + int nmax; + double **ftmp; }; } @@ -60,6 +63,14 @@ E: Pair style is incompatible with KSpace style If a pair style with a long-range Coulombic component is selected, then a kspace style must also be used. +E: Must use 'kspace_modify pressure/scalar no' to obtain per-atom virial with kspace_style MSM + +The kspace scalar pressure option cannot be used to obtain per-atom virial. + +E: Must use 'kspace_modify pressure/scalar no' for rRESPA with kspace_style MSM + +The kspace scalar pressure option cannot (yet) be used with rRESPA. + E: Pair cutoff < Respa interior cutoff One or more pairwise cutoffs are too short to use with the specified diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp index fe8e8ce9bc..726f089877 100644 --- a/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp +++ b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.cpp @@ -54,6 +54,9 @@ PairLJSDKCoulMSM::PairLJSDKCoulMSM(LAMMPS *lmp) : PairLJSDKCoulLong(lmp) void PairLJSDKCoulMSM::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' with Pair style"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h index af79b6e0bf..d9d086f9ab 100644 --- a/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h +++ b/src/USER-CG-CMM/pair_lj_sdk_coul_msm.h @@ -46,3 +46,12 @@ private: #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with Pair style + +The kspace scalar pressure option is not (yet) compatible with at least one of +the defined Pair styles. + +*/ diff --git a/src/USER-OMP/msm_cg_omp.cpp b/src/USER-OMP/msm_cg_omp.cpp index 5120cfb318..23f94349a0 100644 --- a/src/USER-OMP/msm_cg_omp.cpp +++ b/src/USER-OMP/msm_cg_omp.cpp @@ -73,6 +73,10 @@ MSMCGOMP::~MSMCGOMP() void MSMCGOMP::compute(int eflag, int vflag) { + if (scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with kspace_style msm/cg/omp"); + const double * const q = atom->q; const int nlocal = atom->nlocal; int i,j,n; diff --git a/src/USER-OMP/msm_cg_omp.h b/src/USER-OMP/msm_cg_omp.h index 1b6472afaa..1036db134f 100644 --- a/src/USER-OMP/msm_cg_omp.h +++ b/src/USER-OMP/msm_cg_omp.h @@ -56,6 +56,10 @@ Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. +E: Must use 'kspace_modify pressure/scalar no' with kspace_style msm/cg/omp + +The kspace scalar pressure option is not compatible with kspace_style msm/cg/omp. + E: Cannot (yet) use MSM with triclinic box This feature is not yet supported. diff --git a/src/USER-OMP/msm_omp.cpp b/src/USER-OMP/msm_omp.cpp index 04b6ae6e46..61780c87a2 100755 --- a/src/USER-OMP/msm_omp.cpp +++ b/src/USER-OMP/msm_omp.cpp @@ -50,6 +50,9 @@ MSMOMP::MSMOMP(LAMMPS *lmp, int narg, char **arg) : void MSMOMP::compute(int eflag, int vflag) { + if (scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with kspace_style msm/omp"); MSM::compute(eflag,vflag); diff --git a/src/USER-OMP/msm_omp.h b/src/USER-OMP/msm_omp.h index a086b7af92..0dc6776c22 100755 --- a/src/USER-OMP/msm_omp.h +++ b/src/USER-OMP/msm_omp.h @@ -44,3 +44,11 @@ namespace LAMMPS_NS { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with kspace_style msm/omp + +The kspace scalar pressure option is not compatible with kspace_style msm/omp. + +*/ \ No newline at end of file diff --git a/src/USER-OMP/pair_born_coul_msm_omp.cpp b/src/USER-OMP/pair_born_coul_msm_omp.cpp index 950e20e326..5ee4c1e376 100755 --- a/src/USER-OMP/pair_born_coul_msm_omp.cpp +++ b/src/USER-OMP/pair_born_coul_msm_omp.cpp @@ -37,6 +37,10 @@ PairBornCoulMSMOMP::PairBornCoulMSMOMP(LAMMPS *lmp) : void PairBornCoulMSMOMP::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with OMP MSM Pair styles"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-OMP/pair_born_coul_msm_omp.h b/src/USER-OMP/pair_born_coul_msm_omp.h index 317279a5c0..09cd1deacb 100755 --- a/src/USER-OMP/pair_born_coul_msm_omp.h +++ b/src/USER-OMP/pair_born_coul_msm_omp.h @@ -46,3 +46,11 @@ class PairBornCoulMSMOMP : public PairBornCoulMSM, public ThrOMP { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with OMP MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with OMP MSM Pair styles. + +*/ \ No newline at end of file diff --git a/src/USER-OMP/pair_buck_coul_msm_omp.cpp b/src/USER-OMP/pair_buck_coul_msm_omp.cpp index 16dcc5de22..129eca130b 100755 --- a/src/USER-OMP/pair_buck_coul_msm_omp.cpp +++ b/src/USER-OMP/pair_buck_coul_msm_omp.cpp @@ -37,6 +37,10 @@ PairBuckCoulMSMOMP::PairBuckCoulMSMOMP(LAMMPS *lmp) : void PairBuckCoulMSMOMP::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with OMP MSM Pair styles"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-OMP/pair_buck_coul_msm_omp.h b/src/USER-OMP/pair_buck_coul_msm_omp.h index 05f2d0a820..07fd57b654 100755 --- a/src/USER-OMP/pair_buck_coul_msm_omp.h +++ b/src/USER-OMP/pair_buck_coul_msm_omp.h @@ -46,3 +46,11 @@ class PairBuckCoulMSMOMP : public PairBuckCoulMSM, public ThrOMP { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with OMP MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with OMP MSM Pair styles. + +*/ \ No newline at end of file diff --git a/src/USER-OMP/pair_coul_msm_omp.cpp b/src/USER-OMP/pair_coul_msm_omp.cpp index 8d965afa62..1af06208ef 100755 --- a/src/USER-OMP/pair_coul_msm_omp.cpp +++ b/src/USER-OMP/pair_coul_msm_omp.cpp @@ -38,6 +38,10 @@ PairCoulMSMOMP::PairCoulMSMOMP(LAMMPS *lmp) : void PairCoulMSMOMP::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' with " + "OMP MSM Pair styles"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-OMP/pair_coul_msm_omp.h b/src/USER-OMP/pair_coul_msm_omp.h index 23e48f251c..29e4632c66 100755 --- a/src/USER-OMP/pair_coul_msm_omp.h +++ b/src/USER-OMP/pair_coul_msm_omp.h @@ -46,3 +46,11 @@ class PairCoulMSMOMP : public PairCoulMSM, public ThrOMP { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with OMP MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with OMP MSM Pair styles. + +*/ \ No newline at end of file diff --git a/src/USER-OMP/pair_lj_charmm_coul_msm_omp.cpp b/src/USER-OMP/pair_lj_charmm_coul_msm_omp.cpp index 1e23b3a062..f13515283a 100755 --- a/src/USER-OMP/pair_lj_charmm_coul_msm_omp.cpp +++ b/src/USER-OMP/pair_lj_charmm_coul_msm_omp.cpp @@ -38,6 +38,10 @@ PairLJCharmmCoulMSMOMP::PairLJCharmmCoulMSMOMP(LAMMPS *lmp) : void PairLJCharmmCoulMSMOMP::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with OMP MSM Pair styles"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-OMP/pair_lj_charmm_coul_msm_omp.h b/src/USER-OMP/pair_lj_charmm_coul_msm_omp.h index 3760f6768b..907912b138 100755 --- a/src/USER-OMP/pair_lj_charmm_coul_msm_omp.h +++ b/src/USER-OMP/pair_lj_charmm_coul_msm_omp.h @@ -46,3 +46,11 @@ class PairLJCharmmCoulMSMOMP : public PairLJCharmmCoulMSM, public ThrOMP { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with OMP MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with OMP MSM Pair styles. + +*/ \ No newline at end of file diff --git a/src/USER-OMP/pair_lj_cut_coul_msm_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_msm_omp.cpp index de6bf676fa..c355d4aff2 100755 --- a/src/USER-OMP/pair_lj_cut_coul_msm_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_msm_omp.cpp @@ -38,6 +38,10 @@ PairLJCutCoulMSMOMP::PairLJCutCoulMSMOMP(LAMMPS *lmp) : void PairLJCutCoulMSMOMP::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with OMP MSM Pair styles"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-OMP/pair_lj_cut_coul_msm_omp.h b/src/USER-OMP/pair_lj_cut_coul_msm_omp.h index 25cc555750..181c4627b8 100755 --- a/src/USER-OMP/pair_lj_cut_coul_msm_omp.h +++ b/src/USER-OMP/pair_lj_cut_coul_msm_omp.h @@ -46,3 +46,11 @@ class PairLJCutCoulMSMOMP : public PairLJCutCoulMSM, public ThrOMP { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with OMP MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with OMP MSM Pair styles. + +*/ \ No newline at end of file diff --git a/src/USER-OMP/pair_lj_sdk_coul_msm_omp.cpp b/src/USER-OMP/pair_lj_sdk_coul_msm_omp.cpp index b46359e0d9..62d1ae56fd 100644 --- a/src/USER-OMP/pair_lj_sdk_coul_msm_omp.cpp +++ b/src/USER-OMP/pair_lj_sdk_coul_msm_omp.cpp @@ -40,6 +40,10 @@ PairLJSDKCoulMSMOMP::PairLJSDKCoulMSMOMP(LAMMPS *lmp) : void PairLJSDKCoulMSMOMP::compute(int eflag, int vflag) { + if (force->kspace->scalar_pressure_flag) + error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' " + "with OMP MSM Pair styles"); + if (eflag || vflag) { ev_setup(eflag,vflag); } else evflag = vflag_fdotr = 0; diff --git a/src/USER-OMP/pair_lj_sdk_coul_msm_omp.h b/src/USER-OMP/pair_lj_sdk_coul_msm_omp.h index c85f7fb438..9e4a922c39 100644 --- a/src/USER-OMP/pair_lj_sdk_coul_msm_omp.h +++ b/src/USER-OMP/pair_lj_sdk_coul_msm_omp.h @@ -47,3 +47,11 @@ class PairLJSDKCoulMSMOMP : public PairLJSDKCoulMSM, public ThrOMP { #endif #endif + +/* ERROR/WARNING messages: + +E: Must use 'kspace_modify pressure/scalar no' with OMP MSM Pair styles + +The kspace scalar pressure option is not (yet) compatible with OMP MSM Pair styles. + +*/ \ No newline at end of file diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp index 15a47073c4..02b5700f6c 100644 --- a/src/compute_pressure.cpp +++ b/src/compute_pressure.cpp @@ -59,8 +59,8 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) : if (icompute < 0) error->all(FLERR,"Could not find compute pressure temperature ID"); if (modify->compute[icompute]->tempflag == 0) - error->all(FLERR, - "Compute pressure temperature ID does not compute temperature"); + error->all(FLERR,"Compute pressure temperature ID does not " + "compute temperature"); } // process optional args @@ -221,6 +221,10 @@ void ComputePressure::compute_vector() if (update->vflag_global != invoked_vector) error->all(FLERR,"Virial was not tallied on needed timestep"); + if (force->kspace && kspace_virial && force->kspace->scalar_pressure_flag) + error->all(FLERR,"Kspace_modify pressure/scalar no required " + "for components of pressure tensor with kspace_style msm"); + // invoke temperature if it hasn't been already double *ke_tensor; diff --git a/src/compute_pressure.h b/src/compute_pressure.h index 8056a1bfaa..093a21c762 100644 --- a/src/compute_pressure.h +++ b/src/compute_pressure.h @@ -74,6 +74,11 @@ E: Compute pressure temperature ID does not compute temperature The compute ID assigned to a pressure computation must compute temperature. +E: Must use 'kspace_modify pressure/scalar no' to get individual +components of pressure tensor with kspace_style MSM + +Self-explanatory. + E: Virial was not tallied on needed timestep You are using a thermo keyword that requires potentials to diff --git a/src/kspace.cpp b/src/kspace.cpp index 5966ee326e..ae56cfb1ad 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -65,6 +65,7 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) slab_volfactor = 1; suffix_flag = Suffix::NONE; adjust_cutoff_flag = 1; + scalar_pressure_flag = 0; accuracy_absolute = -1.0; accuracy_real_6 = -1.0; @@ -490,6 +491,12 @@ void KSpace::modify_params(int narg, char **arg) if (splittol >= 1.0) error->all(FLERR,"Kspace_modify eigtol must be smaller than one"); iarg += 2; + } else if (strcmp(arg[iarg],"pressure/scalar") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); + if (strcmp(arg[iarg+1],"yes") == 0) scalar_pressure_flag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) scalar_pressure_flag = 0; + else error->all(FLERR,"Illegal kspace_modify command"); + iarg += 2; } else error->all(FLERR,"Illegal kspace_modify command"); } } diff --git a/src/kspace.h b/src/kspace.h index ed5b93ee07..e197848d8f 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -48,6 +48,7 @@ class KSpace : protected Pointers { int neighrequest_flag; // used to avoid obsolete construction of neighbor lists int mixflag; // 1 if geometric mixing rules are enforced for LJ coefficients int slabflag; + int scalar_pressure_flag; // 1 if using MSM fast scalar pressure double slab_volfactor;