Commit JT 081418

- initial commit pppm_spin branch
- copied short_range spin files (src/SPIN)
- copied/renamed Stan's file (from pppm_dipole branch)

src/DIPOLE/pair_lj_cut_dipole_long.cpp

@@ -44,7 +44,7 @@ using namespace MathConst;
 PairLJCutDipoleLong::PairLJCutDipoleLong(LAMMPS *lmp) : Pair(lmp)
 {
   single_enable = 0;
-  ewaldflag = dipoleflag = 1;
+  ewaldflag = pppmflag = dipoleflag = 1;
   respa_enable = 0;
 }
 

src/KSPACE/pppm.h

@@ -41,7 +41,7 @@ class PPPM : public KSpace {
   virtual ~PPPM();
   virtual void init();
   virtual void setup();
-  void setup_grid();
+  virtual void setup_grid();
   virtual void compute(int, int);
   virtual int timing_1d(int, double &);
   virtual int timing_3d(int, double &);
@@ -105,10 +105,10 @@ class PPPM : public KSpace {
   double qdist;                // distance from O site to negative charge
   double alpha;                // geometric factor
 
-  void set_grid_global();
+  virtual void set_grid_global();
   void set_grid_local();
   void adjust_gewald();
-  double newton_raphson_f();
+  virtual double newton_raphson_f();
   double derivf();
   double final_accuracy();
 
@@ -145,7 +145,7 @@ class PPPM : public KSpace {
   void compute_drho1d(const FFT_SCALAR &, const FFT_SCALAR &,
                      const FFT_SCALAR &);
   void compute_rho_coeff();
-  void slabcorr();
+  virtual void slabcorr();
 
   // grid communication
 

src/KSPACE/pppm_spin.h

@@ -0,0 +1,213 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef KSPACE_CLASS
+
+KSpaceStyle(pppm/dipole,PPPMDipole)
+
+#else
+
+#ifndef LMP_PPPM_DIPOLE_H
+#define LMP_PPPM_DIPOLE_H
+
+#include "pppm.h"
+
+namespace LAMMPS_NS {
+
+class PPPMDipole : public PPPM {
+ public:
+  PPPMDipole(class LAMMPS *, int, char **);
+  virtual ~PPPMDipole();
+  void init();
+  void setup();
+  void setup_grid();
+  void compute(int, int);
+  int timing_1d(int, double &);
+  int timing_3d(int, double &);
+  double memory_usage();
+
+ protected:
+  void set_grid_global();
+  double newton_raphson_f();
+
+  void allocate();
+  void allocate_peratom();
+  void deallocate();
+  void deallocate_peratom();
+  void compute_gf_denom();
+
+  void slabcorr();
+
+  // grid communication
+
+  void pack_forward(int, FFT_SCALAR *, int, int *);
+  void unpack_forward(int, FFT_SCALAR *, int, int *);
+  void pack_reverse(int, FFT_SCALAR *, int, int *);
+  void unpack_reverse(int, FFT_SCALAR *, int, int *);
+
+  // dipole
+
+  FFT_SCALAR ***densityx_brick_dipole,***densityy_brick_dipole,***densityz_brick_dipole;
+  FFT_SCALAR ***vdxx_brick_dipole,***vdyy_brick_dipole,***vdzz_brick_dipole;
+  FFT_SCALAR ***vdxy_brick_dipole,***vdxz_brick_dipole,***vdyz_brick_dipole;
+  FFT_SCALAR ***ux_brick_dipole,***uy_brick_dipole,***uz_brick_dipole;
+  FFT_SCALAR ***v0x_brick_dipole,***v1x_brick_dipole,***v2x_brick_dipole;
+  FFT_SCALAR ***v3x_brick_dipole,***v4x_brick_dipole,***v5x_brick_dipole;
+  FFT_SCALAR ***v0y_brick_dipole,***v1y_brick_dipole,***v2y_brick_dipole;
+  FFT_SCALAR ***v3y_brick_dipole,***v4y_brick_dipole,***v5y_brick_dipole;
+  FFT_SCALAR ***v0z_brick_dipole,***v1z_brick_dipole,***v2z_brick_dipole;
+  FFT_SCALAR ***v3z_brick_dipole,***v4z_brick_dipole,***v5z_brick_dipole;
+  FFT_SCALAR *work3,*work4;
+  FFT_SCALAR *densityx_fft_dipole,*densityy_fft_dipole,*densityz_fft_dipole;
+  class GridComm *cg_dipole;
+  class GridComm *cg_peratom_dipole;
+  int only_dipole_flag;
+  double musum,musqsum,mu2;
+  double find_gewald_dipole(double, double, bigint, double, double);
+  double newton_raphson_f_dipole(double, double, bigint, double, double);
+  double derivf_dipole(double, double, bigint, double, double);
+  double compute_df_kspace_dipole();
+  double compute_qopt_dipole();
+  void compute_gf_dipole();
+  void make_rho_dipole();
+  void brick2fft_dipole();
+  void poisson_ik_dipole();
+  void poisson_peratom_dipole();
+  void fieldforce_ik_dipole();
+  void fieldforce_peratom_dipole();
+  double final_accuracy_dipole();
+  void musum_musq();
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Cannot (yet) use charges with Kspace style PPPMDipole
+
+Charge-dipole interactions are not yet implemented in PPPMDipole so this
+feature is not yet supported.
+
+E: Must redefine kspace_style after changing to triclinic box
+
+Self-explanatory.
+
+E: Kspace style requires atom attribute mu
+
+The atom style defined does not have this attribute.
+
+E: Cannot (yet) use kspace_modify diff ad with dipoles
+
+This feature is not yet supported.
+
+E: Cannot (yet) use 'electron' units with dipoles
+
+This feature is not yet supported.
+
+E: Cannot yet use triclinic cells with PPPMDipole
+
+This feature is not yet supported.
+
+E: Cannot yet use TIP4P with PPPMDipole
+
+This feature is not yet supported.
+
+E: Cannot use nonperiodic boundaries with PPPM
+
+For kspace style pppm, all 3 dimensions must have periodic boundaries
+unless you use the kspace_modify command to define a 2d slab with a
+non-periodic z dimension.
+
+E: Incorrect boundaries with slab PPPM
+
+Must have periodic x,y dimensions and non-periodic z dimension to use
+2d slab option with PPPM.
+
+E: PPPM order cannot be < 2 or > than %d
+
+This is a limitation of the PPPM implementation in LAMMPS.
+
+E: KSpace style is incompatible with Pair style
+
+Setting a kspace style requires that a pair style with matching
+long-range dipole components be used.
+
+W: Reducing PPPM order b/c stencil extends beyond nearest neighbor processor
+
+This may lead to a larger grid than desired. See the kspace_modify overlap
+command to prevent changing of the PPPM order.
+
+E: PPPM order < minimum allowed order
+
+The default minimum order is 2. This can be reset by the
+kspace_modify minorder command.
+
+E: PPPM grid stencil extends beyond nearest neighbor processor
+
+This is not allowed if the kspace_modify overlap setting is no.
+
+E: KSpace accuracy must be > 0
+
+The kspace accuracy designated in the input must be greater than zero.
+
+E: Could not compute grid size
+
+The code is unable to compute a grid size consistent with the desired
+accuracy. This error should not occur for typical problems. Please
+send an email to the developers.
+
+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 accuracy.
+
+E: Could not compute g_ewald
+
+The Newton-Raphson solver failed to converge to a good value for
+g_ewald. This error should not occur for typical problems. Please
+send an email to the developers.
+
+E: Non-numeric box dimensions - simulation unstable
+
+The box size has apparently blown up.
+
+E: Out of range atoms - cannot compute PPPM
+
+One or more atoms are attempting to map their charge to a PPPM grid
+point that is not owned by a processor. This is likely for one of two
+reasons, both of them bad. First, it may mean that an atom near the
+boundary of a processor's sub-domain has moved more than 1/2 the
+"neighbor skin distance"_neighbor.html without neighbor lists being
+rebuilt and atoms being migrated to new processors. This also means
+you may be missing pairwise interactions that need to be computed.
+The solution is to change the re-neighboring criteria via the
+"neigh_modify"_neigh_modify command. The safest settings are "delay 0
+every 1 check yes". Second, it may mean that an atom has moved far
+outside a processor's sub-domain or even the entire simulation box.
+This indicates bad physics, e.g. due to highly overlapping atoms, too
+large a timestep, etc.
+
+E: Using kspace solver PPPMDipole on system with no dipoles
+
+Must have non-zero dipoles with PPPMDipole.
+
+E: Must use kspace_modify gewald for system with no dipoles
+
+Self-explanatory.
+
+*/

src/SPIN/pair_spin_long.cpp

@@ -0,0 +1,550 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ------------------------------------------------------------------------
+   Contributing authors: Julien Tranchida (SNL)
+                         Stan Moore (SNL)
+   
+   Please cite the related publication:
+   Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018). 
+   Massively parallel symplectic algorithm for coupled magnetic spin dynamics 
+   and molecular dynamics. arXiv preprint arXiv:1801.10233.
+------------------------------------------------------------------------- */
+
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include "pair_spin_long.h"
+#include "atom.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "fix_nve_spin.h"
+#include "force.h"
+#include "kspace.h"
+#include "math_const.h"
+#include "memory.h"
+#include "modify.h"
+#include "error.h"
+#include "update.h"
+
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+
+#define EWALD_F   1.12837917
+#define EWALD_P   0.3275911
+#define A1        0.254829592
+#define A2       -0.284496736
+#define A3        1.421413741
+#define A4       -1.453152027
+#define A5        1.061405429
+
+/* ---------------------------------------------------------------------- */
+
+PairSpinLong::PairSpinLong(LAMMPS *lmp) : PairSpin(lmp),
+lockfixnvespin(NULL)
+{
+  single_enable = 0;
+  ewaldflag = pppmflag = 1;
+  respa_enable = 0;
+  no_virial_fdotr_compute = 1;
+  lattice_flag = 0;
+
+  hbar = force->hplanck/MY_2PI;		// eV/(rad.THz)
+  mub = 5.78901e-5;                	// in eV/T
+  mu_0 = 1.2566370614e-6;		// in T.m/A
+  mub2mu0 = mub * mub * mu_0;		// in eV
+  mub2mu0hbinv = mub2mu0 / hbar;	// in rad.THz
+
+}
+
+/* ----------------------------------------------------------------------
+   free all arrays
+------------------------------------------------------------------------- */
+
+PairSpinLong::~PairSpinLong()
+{
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(cutsq);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairSpinLong::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,inum,jnum,itype,jtype;  
+  double r,rinv,r2inv,rsq;
+  double grij,expm2,t,erfc;
+  double bij[4];
+  double evdwl,ecoul;
+  double xi[3],rij[3];
+  double spi[4],spj[4],fi[3],fmi[3];
+  double pre1,pre2,pre3;
+  int *ilist,*jlist,*numneigh,**firstneigh;  
+
+  evdwl = ecoul = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  double **x = atom->x;
+  double **f = atom->f;
+  double **fm = atom->fm;
+  double **sp = atom->sp;	
+  int *type = atom->type;  
+  int nlocal = atom->nlocal;  
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  pre1 = 2.0 * g_ewald / MY_PIS;
+  pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
+  pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
+
+  // computation of the exchange interaction
+  // loop over atoms and their neighbors
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xi[0] = x[i][0];
+    xi[1] = x[i][1];
+    xi[2] = x[i][2];
+    jlist = firstneigh[i];
+    jnum = numneigh[i]; 
+    spi[0] = sp[i][0]; 
+    spi[1] = sp[i][1]; 
+    spi[2] = sp[i][2];
+    spi[3] = sp[i][3];
+    itype = type[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      j &= NEIGHMASK;
+      jtype = type[j];
+
+      spj[0] = sp[j][0]; 
+      spj[1] = sp[j][1]; 
+      spj[2] = sp[j][2]; 
+      spj[3] = sp[j][3]; 
+
+      evdwl = 0.0;
+
+      fi[0] = fi[1] = fi[2] = 0.0;
+      fmi[0] = fmi[1] = fmi[2] = 0.0;
+      bij[0] = bij[1] = bij[2] = bij[3] = 0.0;
+     
+      rij[0] = x[j][0] - xi[0];
+      rij[1] = x[j][1] - xi[1];
+      rij[2] = x[j][2] - xi[2];
+      rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
+
+      if (rsq < cutsq[itype][jtype]) {
+        r2inv = 1.0/rsq;
+        rinv = sqrt(r2inv);
+
+        if (rsq < cut_spinsq) {
+          r = sqrt(rsq);
+          grij = g_ewald * r;
+          expm2 = exp(-grij*grij);
+          t = 1.0 / (1.0 + EWALD_P*grij);
+          erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+
+          bij[0] = erfc * rinv;
+          bij[1] = (bij[0] + pre1*expm2) * r2inv;
+          bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
+          bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
+
+	  compute_long(i,j,rij,bij,fmi,spi,spj);
+	  compute_long_mech(i,j,rij,bij,fmi,spi,spj);
+          
+	}
+      }
+
+      // force accumulation
+
+      f[i][0] += fi[0] * mub2mu0;	 
+      f[i][1] += fi[1] * mub2mu0;	  	  
+      f[i][2] += fi[2] * mub2mu0;
+      fm[i][0] += fmi[0] * mub2mu0hbinv;	 
+      fm[i][1] += fmi[1] * mub2mu0hbinv;	  	  
+      fm[i][2] += fmi[2] * mub2mu0hbinv;
+
+      if (newton_pair || j < nlocal) {
+	f[j][0] -= fi[0];	 
+        f[j][1] -= fi[1];	  	  
+        f[j][2] -= fi[2];
+      }
+
+      if (eflag) {
+	if (rsq <= cut_spinsq) {
+	  evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] + 
+	    spi[2]*fmi[2];
+	  evdwl *= hbar;
+	}
+      } else evdwl = 0.0;
+
+
+      if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
+	  evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
+
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   update the pair interaction fmi acting on the spin ii
+------------------------------------------------------------------------- */
+
+void PairSpinLong::compute_single_pair(int ii, double fmi[3])
+{
+  int i,j,jj,jnum,itype,jtype;  
+  double r,rinv,r2inv,rsq;
+  double grij,expm2,t,erfc;
+  double bij[4],xi[3],rij[3],spi[4],spj[4];
+  double pre1,pre2,pre3;
+  int *ilist,*jlist,*numneigh,**firstneigh;  
+
+  double **x = atom->x;
+  double **sp = atom->sp;	
+  int *type = atom->type;  
+
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+
+  pre1 = 2.0 * g_ewald / MY_PIS;
+  pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
+  pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
+
+  // computation of the exchange interaction
+  // loop over neighbors of atom i
+    
+  i = ilist[ii];
+  xi[0] = x[i][0];
+  xi[1] = x[i][1];
+  xi[2] = x[i][2];
+  spi[0] = sp[i][0]; 
+  spi[1] = sp[i][1]; 
+  spi[2] = sp[i][2];
+  spi[3] = sp[i][3];
+  jlist = firstneigh[i];
+  jnum = numneigh[i]; 
+  itype = type[i];
+
+  for (jj = 0; jj < jnum; jj++) {
+    j = jlist[jj];
+    j &= NEIGHMASK;
+    jtype = type[j];
+
+    spj[0] = sp[j][0]; 
+    spj[1] = sp[j][1]; 
+    spj[2] = sp[j][2]; 
+    spj[3] = sp[j][3]; 
+
+    fmi[0] = fmi[1] = fmi[2] = 0.0;
+    bij[0] = bij[1] = bij[2] = bij[3] = 0.0;
+   
+    rij[0] = x[j][0] - xi[0];
+    rij[1] = x[j][1] - xi[1];
+    rij[2] = x[j][2] - xi[2];
+    rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
+
+    if (rsq < cutsq[itype][jtype]) {
+      r2inv = 1.0/rsq;
+      rinv = sqrt(r2inv);
+
+      if (rsq < cut_spinsq) {
+        r = sqrt(rsq);
+        grij = g_ewald * r;
+        expm2 = exp(-grij*grij);
+        t = 1.0 / (1.0 + EWALD_P*grij);
+        erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
+
+        bij[0] = erfc * rinv;
+        bij[1] = (bij[0] + pre1*expm2) * r2inv;
+        bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
+        bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
+
+        compute_long(i,j,rij,bij,fmi,spi,spj);
+        
+      }
+    }
+  }
+
+  // force accumulation
+
+  fmi[0] *= mub2mu0hbinv;	 
+  fmi[1] *= mub2mu0hbinv;	  	  
+  fmi[2] *= mub2mu0hbinv;
+}
+
+/* ----------------------------------------------------------------------
+   compute exchange interaction between spins i and j
+------------------------------------------------------------------------- */
+
+void PairSpinLong::compute_long(int i, int j, double rij[3], 
+    double bij[4], double fmi[3], double spi[4], double spj[4])
+{
+  double sjdotr;
+  double b1,b2,gigj;
+
+  gigj = spi[3] * spj[3];
+  sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2];
+
+  b1 = bij[1];
+  b2 = bij[2];
+
+  fmi[0] += gigj * (b2 * sjdotr *rij[0] - b1 * spj[0]);
+  fmi[1] += gigj * (b2 * sjdotr *rij[1] - b1 * spj[1]);
+  fmi[2] += gigj * (b2 * sjdotr *rij[2] - b1 * spj[2]);
+}
+
+/* ----------------------------------------------------------------------
+   compute the mechanical force due to the exchange interaction between atom i and atom j
+------------------------------------------------------------------------- */
+
+void PairSpinLong::compute_long_mech(int i, int j, double rij[3],
+    double bij[4], double fi[3], double spi[3], double spj[3])
+{
+  double sdots,sidotr,sjdotr,b2,b3;
+  double g1,g2,g1b2_g2b3,gigj;
+
+  gigj = spi[3] * spj[3];
+  sdots = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2];
+  sidotr = spi[0]*rij[0] + spi[1]*rij[1] + spi[2]*rij[2];
+  sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2];
+
+  b2 = bij[2];
+  b3 = bij[3];
+  g1 = sdots;
+  g2 = -sidotr*sjdotr;
+  g1b2_g2b3 = g1*b2 + g2*b3;
+
+  fi[0] += gigj * (rij[0] * g1b2_g2b3 + 
+      b2 * (sjdotr*spi[0] + sidotr*spj[0]));
+  fi[1] += gigj * (rij[1] * g1b2_g2b3 + 
+      b2 * (sjdotr*spi[1] + sidotr*spj[1]));
+  fi[2] += gigj * (rij[2] * g1b2_g2b3 + 
+      b2 * (sjdotr*spi[2] + sidotr*spj[2]));
+}
+
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairSpinLong::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  memory->create(setflag,n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  memory->create(cutsq,n+1,n+1,"pair:cutsq");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairSpinLong::settings(int narg, char **arg)
+{
+  if (narg < 1 || narg > 2)
+    error->all(FLERR,"Incorrect args in pair_style command");
+
+  if (strcmp(update->unit_style,"metal") != 0)
+    error->all(FLERR,"Spin simulations require metal unit style");
+
+  cut_spin = force->numeric(FLERR,arg[0]);
+
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairSpinLong::coeff(int narg, char **arg)
+{
+  if (narg < 4 || narg > 5)
+    error->all(FLERR,"Incorrect args for pair coefficients");
+  if (!allocated) allocate();
+
+  // check if args correct
+
+  if (strcmp(arg[2],"long") != 0)
+    error->all(FLERR,"Incorrect args in pair_style command");
+  if (narg != 3) 
+    error->all(FLERR,"Incorrect args in pair_style command");
+  
+  int ilo,ihi,jlo,jhi;
+  force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
+  force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    for (int j = MAX(jlo,i); j <= jhi; j++) {
+      setflag[i][j] = 1;
+      count++;
+    }
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairSpinLong::init_one(int i, int j)
+{
+  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+  
+  double cut = cut_spin;
+  return cut;
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairSpinLong::init_style()
+{
+  if (!atom->sp_flag)
+    error->all(FLERR,"Pair spin requires atom/spin style");
+  
+  // need a full neighbor list
+
+  int irequest = neighbor->request(this,instance_me);
+  neighbor->requests[irequest]->half = 0;
+  neighbor->requests[irequest]->full = 1;
+  
+  // checking if nve/spin is a listed fix
+
+  int ifix = 0;
+  while (ifix < modify->nfix) {
+    if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
+    ifix++;
+  }
+  if (ifix == modify->nfix)
+    error->all(FLERR,"pair/spin style requires nve/spin");
+
+  // get the lattice_flag from nve/spin
+
+  for (int i = 0; i < modify->nfix; i++) {
+    if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
+      lockfixnvespin = (FixNVESpin *) modify->fix[i];
+      lattice_flag = lockfixnvespin->lattice_flag;
+    }
+  }
+
+  // insure use of KSpace long-range solver, set g_ewald
+
+  if (force->kspace == NULL)
+    error->all(FLERR,"Pair style requires a KSpace style");
+
+  g_ewald = force->kspace->g_ewald;
+
+  cut_spinsq = cut_spin * cut_spin;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairSpinLong::write_restart(FILE *fp)
+{
+  write_restart_settings(fp);
+
+  int i,j;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      fwrite(&setflag[i][j],sizeof(int),1,fp);
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairSpinLong::read_restart(FILE *fp)
+{
+  read_restart_settings(fp);
+
+  allocate();
+
+  int i,j;
+  int me = comm->me;
+  for (i = 1; i <= atom->ntypes; i++)
+    for (j = i; j <= atom->ntypes; j++) {
+      if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
+      MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
+    }
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to restart file
+------------------------------------------------------------------------- */
+
+void PairSpinLong::write_restart_settings(FILE *fp)
+{
+  fwrite(&cut_spin,sizeof(double),1,fp);
+  fwrite(&mix_flag,sizeof(int),1,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads from restart file, bcasts
+------------------------------------------------------------------------- */
+
+void PairSpinLong::read_restart_settings(FILE *fp)
+{
+  if (comm->me == 0) {
+    fread(&cut_spin,sizeof(double),1,fp);
+    fread(&mix_flag,sizeof(int),1,fp);
+  }
+  MPI_Bcast(&cut_spin,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void *PairSpinLong::extract(const char *str, int &dim)
+{
+  if (strcmp(str,"cut") == 0) {
+    dim = 0;
+    return (void *) &cut_spin;
+  } else if (strcmp(str,"cut_coul") == 0) {
+    dim = 0;
+    return (void *) &cut_spin;
+  } else if (strcmp(str,"ewald_order") == 0) {
+    ewald_order = 0;
+    ewald_order |= 1<<1;
+    ewald_order |= 1<<3;
+    dim = 0;
+    return (void *) &ewald_order;
+  } else if (strcmp(str,"ewald_mix") == 0) {
+    dim = 0;
+    return (void *) &mix_flag;
+  }
+  return NULL;
+}

src/SPIN/pair_spin_long.h

@@ -0,0 +1,97 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(spin/long,PairSpinLong)
+
+#else
+
+#ifndef LMP_PAIR_SPIN_LONG_H
+#define LMP_PAIR_SPIN_LONG_H
+
+#include "pair_spin.h"
+
+namespace LAMMPS_NS {
+
+class PairSpinLong : public PairSpin {
+ public:
+  double cut_coul;
+  double **sigma;
+
+  PairSpinLong(class LAMMPS *);
+  ~PairSpinLong();
+  void settings(int, char **);
+  void coeff(int, char **);
+  double init_one(int, int);
+  void init_style();
+  void *extract(const char *, int &); 
+  
+  void compute(int, int);
+  void compute_single_pair(int, double *);
+
+  void compute_long(int, int, double *, double *, double *, 
+      double *, double *);
+  void compute_long_mech(int, int, double *, double *, double *, 
+      double *, double *);
+
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_restart_settings(FILE *);
+  void read_restart_settings(FILE *);
+
+ protected:
+  double hbar;	 		// reduced Planck's constant
+  double mub;			// Bohr's magneton
+  double mu_0;			// vacuum permeability
+  double mub2mu0;		// prefactor for mech force
+  double mub2mu0hbinv;		// prefactor for mag force
+  double cut_spin, cut_spinsq;
+
+  double g_ewald;
+  int ewald_order;
+
+  int lattice_flag;			// flag for mech force computation
+  class FixNVESpin *lockfixnvespin;	// ptr for setups
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Incorrect args in pair_style command
+
+Self-explanatory.
+
+E: Incorrect args for pair coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Pair dipole/long requires atom attributes q, mu, torque
+
+The atom style defined does not have these attributes.
+
+E: Cannot (yet) use 'electron' units with dipoles
+
+This feature is not yet supported.
+
+E: Pair style requires a KSpace style
+
+No kspace style is defined.
+
+*/

src/kspace.cpp

@@ -268,7 +268,7 @@ void KSpace::ev_setup(int eflag, int vflag, int alloc)
    called initially, when particle count changes, when charges are changed
 ------------------------------------------------------------------------- */
 
-void KSpace::qsum_qsq()
+void KSpace::qsum_qsq(int warning_flag)
 {
   const double * const q = atom->q;
   const int nlocal = atom->nlocal;
@@ -285,7 +285,7 @@ void KSpace::qsum_qsq()
   MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world);
   MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
 
-  if ((qsqsum == 0.0) && (comm->me == 0) && warn_nocharge) {
+  if ((qsqsum == 0.0) && (comm->me == 0) && warn_nocharge && warning_flag) {
     error->warning(FLERR,"Using kspace solver on system with no charge");
     warn_nocharge = 0;
   }

src/kspace.h

@@ -108,7 +108,7 @@ class KSpace : protected Pointers {
 
   // public so can be called by commands that change charge
 
-  void qsum_qsq();
+  void qsum_qsq(int warning_flag = 1);
 
   // general child-class methods