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@12950 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-01-19 22:14:11 +00:00
parent 3eab216a55
commit 1eb5be1591
23 changed files with 124 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/GPU/pppm_gpu.cpp
View File

@ -290,13 +290,11 @@ void PPPMGPU::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
// sum energy across procs and add in volume-dependent term

12
src/KSPACE/ewald.cpp
View File

@ -119,7 +119,7 @@ void Ewald::init()
scale = 1.0;
qqrd2e = force->qqrd2e;
qsum_qsq(0);
qsum_qsq();
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -431,13 +431,11 @@ void Ewald::compute(int eflag, int vflag)
if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2];
}
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum global energy across Kspace vevs and add in volume-dependent term

58
src/KSPACE/ewald_disp.cpp
View File

@ -140,13 +140,7 @@ void EwaldDisp::init()
pair->init(); // so B is defined
init_coeffs();
init_coeff_sums();
double qsum, qsqsum, bsbsum;
qsum = qsqsum = bsbsum = 0.0;
if (function[0]) {
qsum = sum[0].x;
qsqsum = sum[0].x2;
}
if (function[0]) qsum_qsq();
// turn off coulombic if no charge
@ -156,6 +150,7 @@ void EwaldDisp::init()
nsums -= 1;
}
double bsbsum = 0.0;
if (function[1]) bsbsum = sum[1].x2;
if (function[2]) bsbsum = sum[2].x2;
@ -246,9 +241,6 @@ void EwaldDisp::setup()
error->all(FLERR,"KSpace accuracy too low");
}
if (qsum_update_flag)
error->all(FLERR,"Kspace style ewald/disp does not support dynamic charges");
bigint natoms = atom->natoms;
double err;
int kxmax = 1;
@ -297,7 +289,8 @@ void EwaldDisp::setup()
compute RMS accuracy for a dimension
------------------------------------------------------------------------- */
double EwaldDisp::rms(int km, double prd, bigint natoms, double q2, double b2, double M2)
double EwaldDisp::rms(int km, double prd, bigint natoms,
double q2, double b2, double M2)
{
double value = 0.0;
@ -520,11 +513,18 @@ void EwaldDisp::init_coeff_sums()
Sum sum_local[EWALD_MAX_NSUMS];
memset(sum_local, 0, EWALD_MAX_NSUMS*sizeof(Sum));
if (function[0]) { // 1/r
double *q = atom->q, *qn = q+atom->nlocal;
for (double *i=q; i<qn; ++i) {
sum_local[0].x += i[0]; sum_local[0].x2 += i[0]*i[0]; }
}
// now perform qsum and qsq via parent qsum_qsq()
sum_local[0].x = 0.0;
sum_local[0].x2 = 0.0;
//if (function[0]) { // 1/r
// double *q = atom->q, *qn = q+atom->nlocal;
// for (double *i=q; i<qn; ++i) {
// sum_local[0].x += i[0]; sum_local[0].x2 += i[0]*i[0]; }
//}
if (function[1]) { // geometric 1/r^6
int *type = atom->type, *ntype = type+atom->nlocal;
for (int *i=type; i<ntype; ++i) {
@ -557,8 +557,8 @@ void EwaldDisp::init_self()
memset(virial_self, 0, EWALD_NFUNCS*sizeof(double));
if (function[0]) { // 1/r
virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*sum[0].x*sum[0].x;
energy_self[0] = sum[0].x2*qscale*g1/MY_PIS-virial_self[0];
virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*qsum*qsum;
energy_self[0] = qsqsum*qscale*g1/MY_PIS-virial_self[0];
}
if (function[1]) { // geometric 1/r^6
virial_self[1] = MY_PI*MY_PIS*g3/(6.0*volume)*sum[1].x*sum[1].x;
@ -597,7 +597,7 @@ void EwaldDisp::init_self_peratom()
double *qi = atom->q, *qn = qi + nlocal;
for (; qi < qn; qi++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) {
double q = *qi;
*vi = cv*q*sum[0].x;
*vi = cv*q*qsum;
*ei = ce*q*q-vi[0];
}
}
@ -676,6 +676,14 @@ void EwaldDisp::compute(int eflag, int vflag)
compute_ek();
compute_force();
//compute_surface(); // assume conducting metal (tinfoil) boundary conditions
// update qsum and qsqsum, if atom count has changed and energy needed
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
compute_energy();
compute_energy_peratom();
compute_virial();
@ -1196,7 +1204,6 @@ void EwaldDisp::compute_virial_dipole()
for (int n = 0; n < 6; n++)
virial[n] += sum[n];
}
}
void EwaldDisp::compute_virial_peratom()
@ -1325,7 +1332,6 @@ void EwaldDisp::compute_virial_peratom()
}
}
/* ----------------------------------------------------------------------
Slab-geometry correction term to dampen inter-slab interactions between
periodically repeating slabs. Yields good approximation to 2D Ewald if
@ -1343,9 +1349,6 @@ void EwaldDisp::compute_slabcorr()
double zprd = domain->zprd;
int nlocal = atom->nlocal;
double qsum = 0.0;
if (function[0]) qsum = sum[0].x;
double dipole = 0.0;
for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2];
@ -1401,7 +1404,8 @@ void EwaldDisp::compute_slabcorr()
double ffact = qscale * (-4.0*MY_PI/volume);
double **f = atom->f;
for (int i = 0; i < nlocal; i++) f[i][2] += ffact * q[i]*(dipole_all - qsum*x[i][2]);
for (int i = 0; i < nlocal; i++)
f[i][2] += ffact * q[i]*(dipole_all - qsum*x[i][2]);
// add on torque corrections
@ -1416,8 +1420,8 @@ void EwaldDisp::compute_slabcorr()
}
/* ----------------------------------------------------------------------
Newton solver used to find g_ewald for LJ systems
------------------------------------------------------------------------- */
Newton solver used to find g_ewald for LJ systems
------------------------------------------------------------------------- */
double EwaldDisp::NewtonSolve(double x, double Rc,
bigint natoms, double vol, double b2)

14
src/KSPACE/msm.cpp
View File

@ -100,6 +100,7 @@ MSM::MSM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
peratom_allocate_flag = 0;
scalar_pressure_flag = 1;
warn_nonneutral = 0;
order = 10;
}
@ -183,7 +184,7 @@ void MSM::init()
scale = 1.0;
qqrd2e = force->qqrd2e;
qsum_qsq(1);
qsum_qsq();
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -507,7 +508,6 @@ void MSM::compute(int eflag, int vflag)
restriction(n);
}
// compute direct interation for top grid level for nonperiodic
// and for second from top grid level for periodic
@ -565,13 +565,11 @@ void MSM::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum global energy across procs and add in self-energy term

11
src/KSPACE/msm_cg.cpp
View File

@ -247,13 +247,11 @@ void MSMCG::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum global energy across procs and add in self-energy term
@ -301,7 +299,6 @@ void MSMCG::compute(int eflag, int vflag)
}
}
}
}
/* ----------------------------------------------------------------------

12
src/KSPACE/pppm.cpp
View File

@ -257,7 +257,7 @@ void PPPM::init()
scale = 1.0;
qqrd2e = force->qqrd2e;
qsum_qsq(0);
qsum_qsq();
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -666,13 +666,11 @@ void PPPM::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum global energy across procs and add in volume-dependent term

10
src/KSPACE/pppm_cg.cpp
View File

@ -207,13 +207,11 @@ void PPPMCG::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum global energy across procs and add in volume-dependent term

16
src/KSPACE/pppm_disp.cpp
View File

@ -293,7 +293,7 @@ void PPPMDisp::init()
qqrd2e = force->qqrd2e;
natoms_original = atom->natoms;
if (function[0]) qsum_qsq(0);
if (function[0]) qsum_qsq();
// if kspace is TIP4P, extract TIP4P params from pair style
// bond/angle are not yet init(), so insure equilibrium request is valid
@ -1114,19 +1114,17 @@ void PPPMDisp::compute(int eflag, int vflag)
fieldforce_none_ik();
if (evflag_atom) cg_peratom_6->forward_comm(this, FORWARD_IK_PERATOM_NONE);
if (evflag_atom)
cg_peratom_6->forward_comm(this, FORWARD_IK_PERATOM_NONE);
}
if (evflag_atom) fieldforce_none_peratom();
}
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum energy across procs and add in volume-dependent term

10
src/KSPACE/pppm_stagger.cpp
View File

@ -201,13 +201,11 @@ void PPPMStagger::compute(int eflag, int vflag)
stagger += 1.0/float(nstagger);
}
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq();
natoms_original = atom->natoms;
}
// sum global energy across procs and add in volume-dependent term

2
src/QEQ/fix_qeq.cpp
View File

@ -107,7 +107,6 @@ FixQEq::FixQEq(LAMMPS *lmp, int narg, char **arg) :
s_hist[i][j] = t_hist[i][j] = 0;
read_file(arg[7]);
}
/* ---------------------------------------------------------------------- */
@ -262,7 +261,6 @@ void FixQEq::reallocate_matrix()
void FixQEq::init_list(int id, NeighList *ptr)
{
list = ptr;
if (force->kspace) force->kspace->qsum_update_flag = 1;
}
/* ---------------------------------------------------------------------- */

6
src/QEQ/fix_qeq_dynamic.cpp
View File

@ -64,7 +64,8 @@ FixQEqDynamic::FixQEqDynamic(LAMMPS *lmp, int narg, char **arg) :
void FixQEqDynamic::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/dynamic requires atom attribute q");
if (!atom->q_flag)
error->all(FLERR,"Fix qeq/dynamic requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/dynamic group has no atoms");
@ -168,8 +169,7 @@ void FixQEqDynamic::pre_force(int vflag)
}
}
if (force->kspace) force->kspace->setup();
if (force->kspace) force->kspace->qsum_qsq();
}
/* ---------------------------------------------------------------------- */

6
src/QEQ/fix_qeq_point.cpp
View File

@ -45,7 +45,8 @@ FixQEqPoint::FixQEqPoint(LAMMPS *lmp, int narg, char **arg) :
void FixQEqPoint::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/point requires atom attribute q");
if (!atom->q_flag)
error->all(FLERR,"Fix qeq/point requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/point group has no atoms");
@ -83,8 +84,7 @@ void FixQEqPoint::pre_force(int vflag)
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
if (force->kspace) force->kspace->setup();
if (force->kspace) force->kspace->qsum_qsq();
}
/* ---------------------------------------------------------------------- */

9
src/QEQ/fix_qeq_shielded.cpp
View File

@ -45,7 +45,8 @@ FixQEqShielded::FixQEqShielded(LAMMPS *lmp, int narg, char **arg) :
void FixQEqShielded::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/shielded requires atom attribute q");
if (!atom->q_flag)
error->all(FLERR,"Fix qeq/shielded requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/shielded group has no atoms");
@ -63,7 +64,8 @@ void FixQEqShielded::init()
int i;
for (i = 1; i <= ntypes; i++) {
if (gamma[i] == 0.0) error->all(FLERR,"Invalid param file for fix qeq/shielded");
if (gamma[i] == 0.0)
error->all(FLERR,"Invalid param file for fix qeq/shielded");
}
if (strstr(update->integrate_style,"respa"))
@ -126,8 +128,7 @@ void FixQEqShielded::pre_force(int vflag)
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
if (force->kspace) force->kspace->setup();
if (force->kspace) force->kspace->qsum_qsq();
}
/* ---------------------------------------------------------------------- */

9
src/QEQ/fix_qeq_slater.cpp
View File

@ -62,7 +62,8 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) :
void FixQEqSlater::init()
{
if (!atom->q_flag) error->all(FLERR,"Fix qeq/slater requires atom attribute q");
if (!atom->q_flag)
error->all(FLERR,"Fix qeq/slater requires atom attribute q");
ngroup = group->count(igroup);
if (ngroup == 0) error->all(FLERR,"Fix qeq/slater group has no atoms");
@ -75,7 +76,8 @@ void FixQEqSlater::init()
int ntypes = atom->ntypes;
for (int i = 1; i <= ntypes; i++) {
if (zeta[i] == 0.0) error->all(FLERR,"Invalid param file for fix qeq/slater");
if (zeta[i] == 0.0)
error->all(FLERR,"Invalid param file for fix qeq/slater");
}
if (strstr(update->integrate_style,"respa"))
@ -101,8 +103,7 @@ void FixQEqSlater::pre_force(int vflag)
matvecs += CG(b_t, t); // CG on t - parallel
calculate_Q();
if (force->kspace) force->kspace->setup();
if (force->kspace) force->kspace->qsum_qsq();
}
/* ---------------------------------------------------------------------- */

2
src/USER-CUDA/pppm_cuda.cpp
View File

@ -691,7 +691,7 @@ void PPPMCuda::compute(int eflag, int vflag)
// extend size of per-atom arrays if necessary
// and force update of device data, if charges can change
if ((cu_atom->update_nmax)||(old_nmax==0)||qsum_update_flag) {
if ((cu_atom->update_nmax)||(old_nmax==0)||qdynamic) {
memory->destroy(part2grid);
nmax = atom->nmax;
memory->create(part2grid,nmax,3,"pppm:part2grid");

10
src/USER-CUDA/pppm_old.cpp
View File

@ -722,13 +722,11 @@ void PPPMOld::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
// sum global energy across procs and add in volume-dependent term

19
src/USER-FEP/compute_fep.cpp
View File

@ -412,15 +412,9 @@ void ComputeFEP::perturb_params()
if (pairflag) force->pair->reinit();
// when perturbing charge and using kspace,
// need to recompute additional params in kspace->setup()
// first backup the state of kspace->qsum_update_flag
// reset KSpace charges if charges have changed
if (chgflag && force->kspace) {
sys_qsum_update_flag = force->kspace->qsum_update_flag;
force->kspace->qsum_update_flag = 1;
force->kspace->setup();
}
if (chgflag && force->kspace) force->kspace->qsum_qsq();
}
@ -461,11 +455,10 @@ void ComputeFEP::restore_params()
}
if (pairflag) force->pair->reinit();
if (chgflag && force->kspace) {
force->kspace->setup();
// restore kspace->qsum_update_flag to original state
force->kspace->qsum_update_flag = sys_qsum_update_flag;
}
// reset KSpace charges if charges have changed
if (chgflag && force->kspace) force->kspace->qsum_qsq();
}

12
src/USER-FEP/fix_adapt_fep.cpp
View File

@ -172,11 +172,6 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
if (adapt[m].which == PAIR)
memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
// when adapting charge and using kspace,
// need to recompute additional params in kspace->setup()
if (chgflag && force->kspace) force->kspace->qsum_update_flag = 1;
id_fix_diam = id_fix_chg = NULL;
}
@ -194,8 +189,6 @@ FixAdaptFEP::~FixAdaptFEP()
}
delete [] adapt;
if (chgflag && force->kspace) force->kspace->qsum_update_flag = 0;
// check nfix in case all fixes have already been deleted
if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
@ -533,10 +526,9 @@ void FixAdaptFEP::change_settings()
if (anypair) force->pair->reinit();
// re-setup KSpace if it exists and adapting charges
// since charges have changed
// reset KSpace charges if charges have changed
if (chgflag && force->kspace) force->kspace->setup();
if (chgflag && force->kspace) force->kspace->qsum_qsq();
}
/* ----------------------------------------------------------------------

10
src/USER-OMP/ewald_omp.cpp
View File

@ -86,14 +86,12 @@ void EwaldOMP::compute(int eflag, int vflag)
MPI_Allreduce(sfacrl,sfacrl_all,kcount,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(sfacim,sfacim_all,kcount,MPI_DOUBLE,MPI_SUM,world);
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
// (n.b. needs to be done outside of the multi-threaded region)
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
// K-space portion of electric field

10
src/USER-OMP/msm_cg_omp.cpp
View File

@ -249,13 +249,11 @@ void MSMCGOMP::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// update qsum and qsqsum, if needed
// update qsum and qsqsum, if atom count has changed and energy needed
if (eflag_global || eflag_atom) {
if (qsum_update_flag || (atom->natoms != natoms_original)) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
// sum global energy across procs and add in self-energy term

12
src/fix_adapt.cpp
View File

@ -175,8 +175,6 @@ FixAdapt::~FixAdapt()
}
delete [] adapt;
if (chgflag && force->kspace) force->kspace->qsum_update_flag = 0;
// check nfix in case all fixes have already been deleted
if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
@ -353,11 +351,6 @@ void FixAdapt::init()
}
}
// when adapting charge and using kspace,
// need to recompute additional params in kspace->setup()
if (chgflag && force->kspace) force->kspace->qsum_update_flag = 1;
// fixes that store initial per-atom values
if (id_fix_diam) {
@ -504,10 +497,9 @@ void FixAdapt::change_settings()
if (anypair) force->pair->reinit();
// re-setup KSpace if it exists and adapting charges
// since charges have changed
// reset KSpace charges if charges have changed
if (chgflag && force->kspace) force->kspace->setup();
if (chgflag && force->kspace) force->kspace->qsum_qsq();
}
/* ----------------------------------------------------------------------

16
src/kspace.cpp
View File

@ -68,8 +68,7 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
suffix_flag = Suffix::NONE;
adjust_cutoff_flag = 1;
scalar_pressure_flag = 0;
qsum_update_flag = 0;
warn_neutral = 1;
warn_nonneutral = 1;
accuracy_absolute = -1.0;
accuracy_real_6 = -1.0;
@ -259,10 +258,10 @@ void KSpace::ev_setup(int eflag, int vflag)
/* ----------------------------------------------------------------------
compute qsum,qsqsum,q2 and give error/warning if not charge neutral
only called initially or when particle count changes
called initially, when particle count changes, when charges are changed
------------------------------------------------------------------------- */
void KSpace::qsum_qsq(int flag)
void KSpace::qsum_qsq()
{
const double * const q = atom->q;
const int nlocal = atom->nlocal;
@ -285,15 +284,14 @@ void KSpace::qsum_qsq(int flag)
q2 = qsqsum * force->qqrd2e;
// not yet sure of the correction needed for non-neutral systems
// so issue warning or error
if (fabs(qsum) > SMALL) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
if (warn_neutral && (comm->me == 0)) {
if (flag) error->all(FLERR,str);
else error->warning(FLERR,str);
}
warn_neutral = 0;
if (!warn_nonneutral) error->all(FLERR,str);
if (warn_nonneutral == 1 && comm->me == 0) error->warning(FLERR,str);
warn_nonneutral = 2;
}
}

11
src/kspace.h
View File

@ -51,9 +51,12 @@ class KSpace : protected Pointers {
// for LJ coefficients
int slabflag;
int scalar_pressure_flag; // 1 if using MSM fast scalar pressure
int qsum_update_flag; // 1 if setup() needs to call qsum_qsq()
double slab_volfactor;
int warn_nonneutral; // 0 = error if non-neutral system
// 1 = warn once if non-neutral system
// 2 = warn, but already warned
int order,order_6,order_allocated;
double accuracy; // accuracy of KSpace solver (force units)
double accuracy_absolute; // user-specifed accuracy in force units
@ -101,6 +104,10 @@ class KSpace : protected Pointers {
void lamda2xvector(double *, double *);
void kspacebbox(double, double *);
// public so can be called by commands that change charge
void qsum_qsq();
// general child-class methods
virtual void init() = 0;
@ -168,7 +175,6 @@ class KSpace : protected Pointers {
int minorder,overlap_allowed;
int adjust_cutoff_flag;
int suffix_flag; // suffix compatibility flag
int warn_neutral; // warn about non-neutral system if 1
bigint natoms_original;
double scale,qqrd2e;
double qsum,qsqsum,q2;
@ -182,7 +188,6 @@ class KSpace : protected Pointers {
int kewaldflag; // 1 if kspace range set for Ewald sum
int kx_ewald,ky_ewald,kz_ewald; // kspace settings for Ewald sum
void qsum_qsq(int);
void pair_check();
void ev_setup(int, int);
double estimate_table_accuracy(double, double);