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

src/kspace.cpp

@@ -26,6 +26,8 @@
 
 using namespace LAMMPS_NS;
 
+#define SMALL 0.00001
+
 /* ---------------------------------------------------------------------- */
 
 KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
@@ -243,6 +245,40 @@ 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
+------------------------------------------------------------------------- */
+
+void KSpace::qsum_qsq(int flag)
+{
+  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;
+  MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+  qsum = tmp;
+  MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+  qsqsum = tmp;
+
+  if (qsqsum == 0.0)
+    error->all(FLERR,"Cannot use kspace solver on system with no charge");
+
+  q2 = qsqsum * force->qqrd2e;
+
+  // 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);
+    if (flag) error->all(FLERR,str);
+    else if (comm->me == 0) error->warning(FLERR,str);
+  }
+}
+
 /* ----------------------------------------------------------------------
    estimate the accuracy of the short-range coulomb tables
 ------------------------------------------------------------------------- */
@@ -335,7 +371,8 @@ void KSpace::lamda2xvector(double *lamda, double *v)
    coords and return the tight (axis-aligned) bounding box, does not 
    preserve vector magnitude
    see http://www.loria.fr/~shornus/ellipsoid-bbox.html and
-   http://yiningkarlli.blogspot.com/2013/02/bounding-boxes-for-ellipsoidsfigure.html
+   http://yiningkarlli.blogspot.com/2013/02/
+     bounding-boxes-for-ellipsoidsfigure.html
 ------------------------------------------------------------------------- */
 
 void KSpace::kspacebbox(double r, double *b)

src/kspace.h

@@ -45,8 +45,10 @@ class KSpace : protected Pointers {
   int tip4pflag;                 // 1 if a TIP4P solver
   int dipoleflag;                // 1 if a dipole solver
   int differentiation_flag;
-  int neighrequest_flag;         // used to avoid obsolete construction of neighbor lists
+  int neighrequest_flag;         // used to avoid obsolete construction 
-  int mixflag;                   // 1 if geometric mixing rules are enforced for LJ coefficients
+                                 // 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;
@@ -57,8 +59,10 @@ class KSpace : protected Pointers {
   double accuracy_absolute;         // user-specifed accuracy in force units
   double accuracy_relative;         // user-specified dimensionless accuracy
                                     // accurary = acc_rel * two_charge_force
-  double accuracy_real_6;           // real space accuracy for dispersion solver (force units)
+  double accuracy_real_6;           // real space accuracy for 
-  double accuracy_kspace_6;         // reciprocal space accuracy for dispersion solver (force units)
+                                    // dispersion solver (force units)
+  double accuracy_kspace_6;         // reciprocal space accuracy for 
+                                    // dispersion solver (force units)
   double two_charge_force;          // force in user units of two point
                                     // charges separated by 1 Angstrom
 
@@ -77,7 +81,7 @@ class KSpace : protected Pointers {
   int collective_flag;            // 1 if use MPI collectives for FFT/remap
   int stagger_flag;               // 1 if using staggered PPPM grids
 
-  double splittol;                // tolerance for when to truncate the splitting
+  double splittol;                // tolerance for when to truncate splitting
 
   KSpace(class LAMMPS *, int, char **);
   virtual ~KSpace();
@@ -116,19 +120,19 @@ class KSpace : protected Pointers {
    see Eq 4 from Parallel Computing 35 (2009) 164–177
 ------------------------------------------------------------------------- */
 
-  double gamma(const double &rho) const {
+  double gamma(const double &rho) const 
+  {
     if (rho <= 1.0) {
       const int split_order = order/2;
       const double rho2 = rho*rho;
       double g = gcons[split_order][0];
       double rho_n = rho2;
-      for (int n=1; n<=split_order; n++) {
+      for (int n = 1; n <= split_order; n++) {
         g += gcons[split_order][n]*rho_n;
         rho_n *= rho2;
       }
       return g;
-    } else
+    } else return (1.0/rho);
-      return (1.0/rho);
   }
 
 /* ----------------------------------------------------------------------
@@ -136,19 +140,19 @@ class KSpace : protected Pointers {
    see Eq 4 from Parallel Computing 35 (2009) 164-177
 ------------------------------------------------------------------------- */
 
-  double dgamma(const double &rho) const {
+  double dgamma(const double &rho) const 
+  {
     if (rho <= 1.0) {
       const int split_order = order/2;
       const double rho2 = rho*rho;
       double dg = dgcons[split_order][0]*rho;
       double rho_n = rho*rho2;
-      for (int n=1; n<split_order; n++) {
+      for (int n = 1; n < split_order; n++) {
         dg += dgcons[split_order][n]*rho_n;
         rho_n *= rho2;
       }
       return dg;
-    } else
+    } else return (-1.0/rho/rho);
-      return (-1.0/rho/rho);
   }
   
   double **get_gcons() { return gcons; }
@@ -160,7 +164,9 @@ class KSpace : protected Pointers {
   int minorder,overlap_allowed;
   int adjust_cutoff_flag;
   int suffix_flag;                  // suffix compatibility flag
-  double scale;
+  bigint natoms_original;
+  double scale,qqrd2e;
+  double qsum,qsqsum,q2;
   double **gcons,**dgcons;          // accumulated per-atom energy/virial
 
   int evflag,evflag_atom;
@@ -171,6 +177,7 @@ 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);