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

src/KSPACE/pppm.cpp

@@ -60,8 +60,8 @@ PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
 {
   if (narg < 1) error->all(FLERR,"Illegal kspace_style pppm command");
 
-  precision = atof(arg[0]);
+  accuracy_relative = atof(arg[0]);
-  
+
   nfactors = 3;
   factors = new int[nfactors];
   factors[0] = 2;
@@ -222,6 +222,11 @@ void PPPM::init()
     error->warning(FLERR,str);
   }
 
+  // set accuracy (force units) from accuracy_relative or accuracy_absolute
+  
+  if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute;
+  else accuracy = accuracy_relative * two_charge_force;
+  
   // setup FFT grid resolution and g_ewald
   // normally one iteration thru while loop is all that is required
   // if grid stencil extends beyond neighbor proc, reduce order and try again
@@ -229,7 +234,6 @@ void PPPM::init()
   int iteration = 0;
 
   while (order > 0) {
-
     if (iteration && me == 0)
       error->warning(FLERR,"Reducing PPPM order b/c stencil extends "
 		     "beyond neighbor processor");
@@ -970,9 +974,7 @@ void PPPM::set_grid()
   acons[7][5] = 1755948832039.0 / 36229939200000.0;
   acons[7][6] = 4887769399.0 / 37838389248.0;
 
-  //double q2 = qsqsum * force->qqrd2e / force->dielectric;
+  double q2 = qsqsum * force->qqrd2e / force->dielectric;
-  double q2 = qsqsum / force->dielectric;
-  bigint natoms = atom->natoms;
 
   // use xprd,yprd,zprd even if triclinic so grid size is the same
   // adjust z dimension for 2d slab PPPM
@@ -984,20 +986,21 @@ void PPPM::set_grid()
   double zprd_slab = zprd*slab_volfactor;
   
   // make initial g_ewald estimate
-  // based on desired error and real space cutoff
+  // based on desired accuracy and real space cutoff
   // fluid-occupied volume used to estimate real-space error
   // zprd used rather than zprd_slab
 
   double h_x,h_y,h_z;
+  bigint natoms = atom->natoms;
 
   if (!gewaldflag)
-    g_ewald = sqrt(-log(precision*sqrt(natoms*cutoff*xprd*yprd*zprd) / 
+    g_ewald = sqrt(-log(accuracy*sqrt(natoms*cutoff*xprd*yprd*zprd) / 
 			(2.0*q2))) / cutoff;
 
-  // set optimal nx_pppm,ny_pppm,nz_pppm based on order and precision
+  // set optimal nx_pppm,ny_pppm,nz_pppm based on order and accuracy
   // nz_pppm uses extended zprd_slab instead of zprd
   // h = 1/g_ewald is upper bound on h such that h*g_ewald <= 1
-  // reduce it until precision target is met
+  // reduce it until accuracy target is met
 
   if (!gridflag) {
     double err;
@@ -1008,21 +1011,21 @@ void PPPM::set_grid()
     nz_pppm = static_cast<int> (zprd_slab/h_z + 1);
 
     err = rms(h_x,xprd,natoms,q2,acons);
-    while (err > precision) {
+    while (err > accuracy) {
       err = rms(h_x,xprd,natoms,q2,acons);
       nx_pppm++;
       h_x = xprd/nx_pppm;
     }
 
     err = rms(h_y,yprd,natoms,q2,acons);
-    while (err > precision) {
+    while (err > accuracy) {
       err = rms(h_y,yprd,natoms,q2,acons);
       ny_pppm++;
       h_y = yprd/ny_pppm;
     }
 
     err = rms(h_z,zprd_slab,natoms,q2,acons);
-    while (err > precision) {
+    while (err > accuracy) {
       err = rms(h_z,zprd_slab,natoms,q2,acons);
       nz_pppm++;
       h_z = zprd_slab/nz_pppm;
@@ -1067,7 +1070,7 @@ void PPPM::set_grid()
     }
   }
 
-  // final RMS precision
+  // final RMS accuracy
 
   double lprx = rms(h_x,xprd,natoms,q2,acons);
   double lpry = rms(h_y,yprd,natoms,q2,acons);
@@ -1092,14 +1095,18 @@ void PPPM::set_grid()
       fprintf(screen,"  G vector (1/distance)= %g\n",g_ewald);
       fprintf(screen,"  grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
       fprintf(screen,"  stencil order = %d\n",order);
-      fprintf(screen,"  RMS precision = %g\n",MAX(lpr,spr));
+      fprintf(screen,"  absolute RMS force accuracy = %g\n",MAX(lpr,spr));
+      fprintf(screen,"  relative force accuracy = %g\n",
+	      MAX(lpr,spr)/two_charge_force);
       fprintf(screen,"  using %s precision FFTs\n",fft_prec);
     }
     if (logfile) {
       fprintf(logfile,"  G vector (1/distance) = %g\n",g_ewald);
       fprintf(logfile,"  grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
       fprintf(logfile,"  stencil order = %d\n",order);
-      fprintf(logfile,"  RMS precision = %g\n",MAX(lpr,spr));
+      fprintf(logfile,"  absolute RMS force accuracy = %g\n",MAX(lpr,spr));
+      fprintf(logfile,"  relative force accuracy = %g\n",
+	      MAX(lpr,spr)/two_charge_force);
       fprintf(logfile,"  using %s precision FFTs\n",fft_prec);
     }
   }
@@ -1128,7 +1135,7 @@ int PPPM::factorable(int n)
 }
 
 /* ----------------------------------------------------------------------
-   compute RMS precision for a dimension
+   compute RMS accuracy for a dimension
 ------------------------------------------------------------------------- */
 
 double PPPM::rms(double h, double prd, bigint natoms,
@@ -1143,7 +1150,7 @@ double PPPM::rms(double h, double prd, bigint natoms,
 }
 
 /* ----------------------------------------------------------------------
-   compute difference in real-space and kspace RMS precision
+   compute difference in real-space and KSpace RMS accuracy
 ------------------------------------------------------------------------- */
 
 double PPPM::diffpr(double h_x, double h_y, double h_z, double q2, 

src/KSPACE/pppm.h

@@ -47,7 +47,6 @@ class PPPM : public KSpace {
 
  protected:
   int me,nprocs;
-  double precision;
   int nfactors;
   int *factors;
   double qsum,qsqsum;
@@ -221,7 +220,7 @@ E: PPPM grid is too large
 
 The global PPPM grid is larger than OFFSET in one or more dimensions.
 OFFSET is currently set to 4096.  You likely need to decrease the
-requested precision.
+requested accuracy.
 
 E: PPPM order has been reduced to 0
 

src/USER-EWALDN/ewald_n.h

@@ -1,82 +1,91 @@
-/* ----------------------------------------------------------------------
+/* ----------------------------------------------------------------------
-   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
-   http://lammps.sandia.gov, Sandia National Laboratories
+   http://lammps.sandia.gov, Sandia National Laboratories
-   Steve Plimpton, sjplimp@sandia.gov
+   Steve Plimpton, sjplimp@sandia.gov
-
+
-   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
-   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
-   certain rights in this software.  This software is distributed under 
+   certain rights in this software.  This software is distributed under 
-   the GNU General Public License.
+   the GNU General Public License.
-
+
-   See the README file in the top-level LAMMPS directory.
+   See the README file in the top-level LAMMPS directory.
-------------------------------------------------------------------------- */
+------------------------------------------------------------------------- */
-
+
-#ifdef KSPACE_CLASS
+#ifdef KSPACE_CLASS
-
+
-KSpaceStyle(ewald/n,EwaldN)
+KSpaceStyle(ewald/n,EwaldN)
-
+
-#else
+#else
-
+
-#ifndef LMP_EWALD_N_H
+#ifndef LMP_EWALD_N_H
-#define LMP_EWALD_N_H
+#define LMP_EWALD_N_H
-
+
-#include "kspace.h"
+#include "kspace.h"
-#include "math_complex.h"
+#include "math_complex.h"
-
+
-namespace LAMMPS_NS {
+namespace LAMMPS_NS {
-
+
-#define EWALD_NORDER	6
+#define EWALD_NORDER	6
-#define EWALD_NFUNCS	4
+#define EWALD_NFUNCS	4
-#define EWALD_MAX_NSUMS	10
+#define EWALD_MAX_NSUMS	10
-#define EWALD_NSUMS	{1, 1, 7, 1}
+#define EWALD_NSUMS	{1, 1, 7, 1}
-
+
-typedef struct cvector { complex x, y, z; } cvector;
+typedef struct cvector { complex x, y, z; } cvector;
-typedef struct hvector { double x, y, z; } hvector;
+typedef struct hvector { double x, y, z; } hvector;
-typedef struct kvector { long x, y, z; } kvector;
+typedef struct kvector { long x, y, z; } kvector;
-
+
-class EwaldN : public KSpace {
+class EwaldN : public KSpace {
- public:
+ public:
-  EwaldN(class LAMMPS *, int, char **);
+  EwaldN(class LAMMPS *, int, char **);
-  ~EwaldN();
+  ~EwaldN();
-  void init();
+  void init();
-  void setup();
+  void setup();
-  void compute(int, int);
+  void compute(int, int);
-  double memory_usage() {return bytes;}
+  double memory_usage() {return bytes;}
-
+
- private:
+ private:
-  double unit[6];
+  double unit[6];
-  int function[EWALD_NFUNCS], first_output;
+  int function[EWALD_NFUNCS], first_output;
-
+
-  int nkvec, nkvec_max, nevec, nevec_max,
+  int nkvec, nkvec_max, nevec, nevec_max,
-      nbox, nfunctions, nsums, sums;
+      nbox, nfunctions, nsums, sums;
-  double bytes;
+  int peratom_allocate_flag;
-  double precision, g2_max;
+  double bytes;
-  double *kenergy, energy_self[EWALD_NFUNCS];
+  double precision, g2_max;
-  double *kvirial, virial_self[EWALD_NFUNCS];
+  double *kenergy, energy_self[EWALD_NFUNCS];
-  cvector *ekr_local;
+  double *kvirial, virial_self[EWALD_NFUNCS];
-  hvector *hvec;
+  double **energy_self_peratom;
-  kvector *kvec;
+  double **virial_self_peratom;
-
+  cvector *ekr_local;
-  double mumurd2e, dielectric, *B, volume;
+  hvector *hvec;
-  struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS];
+  kvector *kvec;
-  complex *cek_local, *cek_global;
+
- 
+  double mumurd2e, dielectric, *B, volume;
-  void reallocate();
+  struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS];
-  void reallocate_atoms();
+  complex *cek_local, *cek_global;
-  void deallocate();
+ 
-  void coefficients();
+  void reallocate();
-  void init_coeffs();
+  void allocate_peratom();
-  void init_coeff_sums();
+  void reallocate_atoms();
-  void init_self();
+  void deallocate();
-  void compute_ek();
+  void deallocate_peratom();
-  void compute_force();
+  void coefficients();
-  void compute_surface();
+  void init_coeffs();
-  void compute_energy(int);
+  void init_coeff_sums();
-  void compute_virial(int);
+  void init_self();
-  void compute_slabcorr(int);
+  void init_self_peratom();
-};
+  void compute_ek();
-
+  void compute_force();
-}
+  void compute_surface();
-
+  void compute_surface_peratom();
-#endif
+  void compute_energy();
-#endif
+  void compute_energy_peratom();
+  void compute_virial();
+  void compute_virial_peratom();
+  void compute_slabcorr();
+};
+
+}
+
+#endif
+#endif

src/fix_box_relax.cpp

@@ -324,7 +324,8 @@ void FixBoxRelax::init()
   temperature = modify->compute[icompute];
 
   icompute = modify->find_compute(id_press);
-  if (icompute < 0) error->all(FLERR,"Pressure ID for fix box/relax does not exist");
+  if (icompute < 0) 
+    error->all(FLERR,"Pressure ID for fix box/relax does not exist");
   pressure = modify->compute[icompute];
 
   pv2e = 1.0 / force->nktv2p;

src/kspace.cpp

@@ -16,6 +16,7 @@
 #include "kspace.h"
 #include "atom.h"
 #include "comm.h"
+#include "force.h"
 #include "memory.h"
 #include "error.h"
 
@@ -32,6 +33,11 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   slabflag = 0;
   slab_volfactor = 1;
 
+  accuracy_absolute = -1.0;
+  two_charge_force = force->qqr2e * 
+    (force->qelectron * force->qelectron) / 
+    (force->angstrom * force->angstrom);
+
   maxeatom = maxvatom = 0;
   eatom = NULL;
   vatom = NULL;
@@ -121,6 +127,10 @@ void KSpace::modify_params(int narg, char **arg)
       if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
       order = atoi(arg[iarg+1]);
       iarg += 2;
+    } else if (strcmp(arg[iarg],"force") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
+      accuracy_absolute = atof(arg[iarg+1]);
+      iarg += 2;
     } else if (strcmp(arg[iarg],"gewald") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command");
       g_ewald = atof(arg[iarg+1]);

src/kspace.h

@@ -45,6 +45,13 @@ class KSpace : protected Pointers {
   int slabflag;
   double scale;
   double slab_volfactor;
+ 
+  double accuracy;                  // accuracy of KSpace solver (force units)
+  double accuracy_absolute;         // user-specifed accuracy in force units
+  double accuracy_relative;         // user-specified dimensionless accuracy
+                                    // accurary = acc_rel * two_charge_force
+  double two_charge_force;          // force in user units of two point
+                                    // charges separated by 1 Angstrom
 
   int evflag,evflag_atom;
   int eflag_either,eflag_global,eflag_atom;