convert to KSspace style rather than fix

examples/scafacos/in.scafacos

@@ -5,17 +5,32 @@ atom_style  full
 
 read_data data.NaCl
 
-velocity	all set 0.0 0.0 0.0 mom no
+replicate 8 8 8
 
-pair_style	zero 1.0
-pair_coeff	* * 
+velocity	all create 1.5 49893
 
 neighbor	1.0 bin
 neigh_modify	delay 0
 
-# using the fmm solver with a tolerance in energy of 10^-3
-fix		1 all scafacos fmm tolerance energy 0.001
+fix             1 all nve
+
+# LAMMPS computes pairwise and long-range Coulombics
+
+#pair_style      coul/long 3.0
+#pair_coeff      * *
+#kspace_style    pppm 1.0e-3
+
+# Scafacos computes entire long-range Coulombics
+# use dummy pair style to perform atom sorting
+
+pair_style	zero 1.0
+pair_coeff	* * 
+
+#fix		2 all scafacos p3m tolerance field 0.001
+
+kspace_style    scafacos p3m tolerance field 0.001
 
 timestep	0.005
+thermo          10
 
 run		100

lib/scafacos/Install.py

+#!/usr/bin/env python
+
+# Install.py tool to download, unpack, build, and link to the Scafacos library
+# used to automate the steps described in the README file in this dir
+
+from __future__ import print_function
+import sys,os,re,subprocess
+
+# help message
+
+help = """
+Syntax from src dir: make lib-scafacos args="-b"
+                 or: make lib-scafacos args="-p /usr/local/scafacos"
+Syntax from lib dir: python Install.py -b
+                 or: python Install.py -p /usr/local/scafacos
+
+specify zero or more options, order does not matter
+
+  -b = download and build the Scafacos library
+  -p = specify folder of existing Scafacos installation
+
+   always creates includelink, liblink to Scafacos dirs
+
+Example:
+
+make lib-scafacos args="-b"   # download/build in lib/scafacos/scafacos
+make lib-scafacos args="-p $HOME/scafacos" # use existing Scafacos installation in $HOME
+"""
+
+# settings
+
+#version = "voro++-0.4.6"
+#url = "http://math.lbl.gov/voro++/download/dir/%s.tar.gz" % version
+
+# print error message or help
+
+def error(str=None):
+  if not str: print(help)
+  else: print("ERROR",str)
+  sys.exit()
+
+# expand to full path name
+# process leading '~' or relative path
+
+def fullpath(path):
+  return os.path.abspath(os.path.expanduser(path))
+
+def which(program):
+  def is_exe(fpath):
+    return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
+
+  fpath, fname = os.path.split(program)
+  if fpath:
+    if is_exe(program):
+      return program
+  else:
+    for path in os.environ["PATH"].split(os.pathsep):
+      path = path.strip('"')
+      exe_file = os.path.join(path, program)
+      if is_exe(exe_file):
+        return exe_file
+
+  return None
+
+def geturl(url,fname):
+  success = False
+
+  if which('curl') != None:
+    cmd = 'curl -L -o "%s" %s' % (fname,url)
+    try:
+      subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+      success = True
+    except subprocess.CalledProcessError as e:
+      print("Calling curl failed with: %s" % e.output.decode('UTF-8'))
+
+  if not success and which('wget') != None:
+    cmd = 'wget -O "%s" %s' % (fname,url)
+    try:
+      subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+      success = True
+    except subprocess.CalledProcessError as e:
+      print("Calling wget failed with: %s" % e.output.decode('UTF-8'))
+
+  if not success:
+    error("Failed to download source code with 'curl' or 'wget'")
+  return
+
+# parse args
+
+args = sys.argv[1:]
+nargs = len(args)
+#if nargs == 0: error()
+
+homepath = "."
+
+buildflag = False
+pathflag = False
+linkflag = True
+
+iarg = 0
+while iarg < nargs:
+  if args[iarg] == "-v":
+    if iarg+2 > nargs: error()
+    version = args[iarg+1]
+    iarg += 2
+  elif args[iarg] == "-p":
+    if iarg+2 > nargs: error()
+    scafacospath = fullpath(args[iarg+1])
+    pathflag = True
+    iarg += 2
+  elif args[iarg] == "-b":
+    buildflag = True
+    iarg += 1
+  else: error()
+
+homepath = fullpath(homepath)
+#homedir = "%s/%s" % (homepath,version)
+homedir = homepath
+
+if (pathflag):
+    if not os.path.isdir(scafacospath): error("Scafacos path does not exist")
+    homedir =scafacospath
+
+if (buildflag and pathflag):
+    error("Cannot use -b and -p flag at the same time")
+
+#if (not buildflag and not pathflag):
+#    error("Have to use either -b or -p flag")
+
+# download and unpack Scafacos tarball
+
+if buildflag:
+  print("Downloading Scafacos ...")
+  geturl(url,"%s/%s.tar.gz" % (homepath,version))
+
+  print("Unpacking Voro++ tarball ...")
+  if os.path.exists("%s/%s" % (homepath,version)):
+    cmd = 'rm -rf "%s/%s"' % (homepath,version)
+    subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+  cmd = 'cd "%s"; tar -xzvf %s.tar.gz' % (homepath,version)
+  subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+  os.remove("%s/%s.tar.gz" % (homepath,version))
+  if os.path.basename(homedir) != version:
+    if os.path.exists(homedir):
+      cmd = 'rm -rf "%s"' % homedir
+      subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+    os.rename("%s/%s" % (homepath,version),homedir)
+
+# build Scafacos
+
+if buildflag:
+  print("Building Scafacos ...")
+  cmd = 'cd "%s"; make CXX=g++ CFLAGS="-fPIC -O3"' % homedir
+  txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+  print(txt.decode('UTF-8'))
+
+# create 2 links in lib/scafacos to Scafacos include/lib dirs
+
+if linkflag:
+  print("Creating links to Scafacos include and lib files")
+  if os.path.isfile("includelink") or os.path.islink("includelink"):
+    os.remove("includelink")
+  if os.path.isfile("liblink") or os.path.islink("liblink"):
+    os.remove("liblink")
+  cmd = 'ln -s "scafacos/include" includelink'
+  subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
+  cmd = 'ln -s "scafacos/lib" liblink'
+  subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)

src/Makefile

@@ -62,8 +62,8 @@ PACKUSER = user-atc user-awpmd user-bocs user-cgdna user-cgsdk user-colvars \
 	   user-intel user-lb user-manifold user-meamc user-meso \
 	   user-mgpt user-misc user-mofff user-molfile \
 	   user-netcdf user-omp user-phonon user-qmmm user-qtb \
-	   user-quip user-reaxc user-smd user-smtbq user-sph user-tally \
-	   user-uef user-vtk
+	   user-quip user-reaxc user-scafacos user-smd user-smtbq \
+	   user-sph user-tally user-uef user-vtk
 
 PACKLIB = compress gpu kim kokkos latte meam mpiio mscg poems \
 	  python reax voronoi \

src/USER-SCAFACOS/Install.sh

@@ -37,31 +37,18 @@ done
 if (test $1 = 1) then
 
   if (test -e ../Makefile.package) then
+    sed -i -e 's/`.*scafacos.*` //' ../Makefile.package
     sed -i -e 's/[^ \t]*scafacos[^ \t]* //' ../Makefile.package
     sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/scafacos/includelink |' ../Makefile.package
     sed -i -e 's|^PKG_PATH =[ \t]*|&-L../../lib/scafacos/liblink |' ../Makefile.package
-    #sed -i -e 's|^PKG_LIB =[ \t]*|&-lscafacos |' ../Makefile.package
-    sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(scafacos_SYSINC) |' ../Makefile.package
-    sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(scafacos_SYSLIB) |' ../Makefile.package
-    sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(scafacos_SYSPATH) |' ../Makefile.package
-  fi
-
-  if (test -e ../Makefile.package.settings) then
-    sed -i -e '/^include.*scafacos.*$/d' ../Makefile.package.settings
-    # multiline form needed for BSD sed on Macs
-    sed -i -e '4 i \
-include ..\/..\/lib\/scafacos\/Makefile.lammps
-' ../Makefile.package.settings
+    sed -i -e 's%^PKG_LIB =[ \t]*%&`grep Libs: ../../lib/scafacos/liblink/pkgconfig/scafacos.pc | cut -d " " -f 2-` %' ../Makefile.package
   fi
 
 elif (test $1 = 0) then
 
-  if (test -e ../Makefile.package) then
+  if (test -e ../Makefile.package) then 
+    sed -i -e 's/`.*scafacos.*` //' ../Makefile.package
     sed -i -e 's/[^ \t]*scafacos[^ \t]* //' ../Makefile.package
   fi
 
-  if (test -e ../Makefile.package.settings) then
-    sed -i -e '/^include.*scafacos.*$/d' ../Makefile.package.settings
-  fi
-
 fi

src/USER-SCAFACOS/fix_scafacos.cpp

@@ -17,6 +17,7 @@
 
 #include <cstdio>
 #include <cstring>
+#include <cstdlib>
 #include "fix_scafacos.h"
 #include "atom.h"
 #include "comm.h"
@@ -81,7 +82,7 @@ FixScafacos::FixScafacos(LAMMPS *lmp, int narg, char **arg) :
       else
         error->all(FLERR,"Illegal fix scafacos command");
       ++arg_index;
-        tolerance_value = atof(arg[arg_index]);
+      tolerance_value = atof(arg[arg_index]);
       ++arg_index; 
     }
     else
@@ -439,6 +440,7 @@ bool FixScafacos::check_result(FCSResult result, int comm_rank)
 {
   if (result) 
   {
+    printf("ScaFaCoS Error: Caught error on task %d.\n", comm_rank);
     std::string err_msg;
     std::stringstream ss;
 
@@ -446,8 +448,9 @@ bool FixScafacos::check_result(FCSResult result, int comm_rank)
        << fcs_result_get_message(result) << "\n";
     err_msg = ss.str();
 
-    error->one(FLERR, err_msg.c_str());
+    error -> all(FLERR, err_msg.c_str());
     fcs_result_destroy(result);
   }
   return true;
 }
+ 

src/USER-SCAFACOS/fix_scafacos.h

@@ -111,4 +111,54 @@ class FixScafacos : public Fix {
 
 /* ERROR/WARNING messages:
 
+E: Must use units metal with fix latte command
+
+UNDOCUMENTED
+
+E: Fix latte currently runs only in serial
+
+UNDOCUMENTED
+
+E: LAMMPS is linked against incompatible LATTE library
+
+UNDOCUMENTED
+
+E: Illegal ... command
+
+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: Fix latte does not yet support a LAMMPS calculation of a Coulomb potential
+
+UNDOCUMENTED
+
+E: Could not find fix latte compute ID
+
+UNDOCUMENTED
+
+E: Fix latte compute ID does not compute pe/atom
+
+UNDOCUMENTED
+
+E: Fix latte requires 3d problem
+
+UNDOCUMENTED
+
+E: Fix latte cannot compute Coulomb potential
+
+UNDOCUMENTED
+
+E: Fix latte requires 3d simulation
+
+UNDOCUMENTED
+
+W: Fix latte should come after all other integration fixes
+
+UNDOCUMENTED
+
+E: Internal LATTE problem
+
+UNDOCUMENTED
+
 */

src/USER-SCAFACOS/scafacos.cpp

+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Rene Halver (JSC)
+------------------------------------------------------------------------- */
+
+#include <cstdio>
+#include <cstring>
+#include <cstdlib>
+#include "scafacos.h"
+#include "atom.h"
+#include "comm.h"
+#include "domain.h"
+#include "force.h"
+#include "memory.h"
+#include "error.h"
+
+// ScaFaCoS library
+
+#include <string>
+#include <sstream>
+#include "fcs.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+Scafacos::Scafacos(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
+{
+  if (narg < 2) error->all(FLERR,"Illegal scafacos command");
+
+  int n = strlen(arg[0]) + 1;
+  method = new char[n];
+  strcpy(method,arg[0]);
+
+  // additional args
+
+  int tflag = 0;
+
+  int iarg = 1;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"tolerance") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal scafacos command");
+      tflag = 1;
+      if (strcmp(arg[iarg+1],"energy") == 0)
+        tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;     
+      else if (strcmp(arg[iarg+1],"energy_rel") == 0)
+        tolerance_type = FCS_TOLERANCE_TYPE_ENERGY_REL;     
+      else if (strcmp(arg[iarg+1],"field") == 0)
+        tolerance_type = FCS_TOLERANCE_TYPE_FIELD;     
+      else if (strcmp(arg[iarg+1],"field_rel") == 0)
+        tolerance_type = FCS_TOLERANCE_TYPE_FIELD_REL;     
+      else if (strcmp(arg[iarg+1],"potential") == 0)
+        tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL;     
+      else if (strcmp(arg[iarg+1],"potential_rel") == 0)
+        tolerance_type = FCS_TOLERANCE_TYPE_POTENTIAL_REL;     
+      else error->all(FLERR,"Illegal scafacos command");
+      tolerance = force->numeric(FLERR,arg[iarg+2]);
+      iarg += 3;
+    } else error->all(FLERR,"Illegal scafacos command");
+  }
+
+  if (!tflag) error->all(FLERR,"Must set tolerance in scafacos command");
+
+  // initializations
+
+  me = comm->me;
+  initialized = 0;
+
+  maxatom = 0;
+  epot = NULL;
+  efield = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+Scafacos::~Scafacos()
+{
+  delete [] method;
+
+  memory->destroy(epot);
+  memory->destroy(efield);
+
+  // NOTE: any clean-up call to ScaFaCoS needed?
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Scafacos::init()
+{
+  // error checks
+  // NOTE: allow triclinic at some point
+
+  if (domain->dimension == 2)
+    error->all(FLERR,"Cannot use ScaFaCoS with 2d simulation");
+
+  if (domain->triclinic)
+    error->all(FLERR,"Cannot use ScaFaCoS with triclinic domain yet");
+
+  // one-time initialization of ScaFaCoS
+
+  if (!initialized) {
+    result = fcs_init(&fcs,method,world);
+    check_result(result);
+
+    setup_handle();
+
+    result = fcs_set_tolerance(fcs,tolerance_type,tolerance);
+    check_result(result);
+
+    if (me == 0) fcs_print_parameters(fcs);
+
+    double **x = atom->x;
+    double *q = atom->q; 
+    int nlocal = atom->nlocal;
+
+    result = fcs_tune(fcs,nlocal,&x[0][0],q);
+  }
+
+  initialized = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void Scafacos::setup() {}
+
+/* ---------------------------------------------------------------------- */
+
+void Scafacos::compute(int eflag, int vflag)
+{
+  double **x = atom->x;
+  double *q = atom->q;
+  int nlocal = atom->nlocal;
+
+  // if simluation box has changed, call fcs_tune()
+
+  if (box_has_changed()) {
+    printf("AAA\n");
+    setup_handle();
+    result = fcs_tune(fcs,nlocal,&x[0][0],q);
+    check_result(result);
+  }
+  
+  // grow epot & efield if necessary
+
+  if (nlocal > maxatom) {
+    memory->destroy(epot);
+    memory->destroy(efield);
+    maxatom = atom->nmax;
+    memory->create(epot,maxatom,"scafacos:epot");
+    memory->create(efield,maxatom,3,"scafacos:efield");
+     
+  }
+
+  // initialize epot & efield
+  // NOTE: is this necessary?
+
+  for (int i = 0; i < nlocal; i++) {
+    epot[i] = 0.0;
+    efield[i][0] = efield[i][1] = efield[i][2] = 0.0;
+  }
+ 
+  // ScaFaCoS calculation of full Coulombics
+
+  result = fcs_run(fcs,nlocal,&x[0][0],q,&efield[0][0],epot);
+  check_result(result);
+
+  // apply Efield to each particle
+  // accumulate total energy
+
+  double **f = atom->f;
+
+  double qone;
+  double myeng = 0.0;
+
+  for (int i = 0; i < nlocal; i++) {
+    qone = q[i];
+    f[i][0] += qone * efield[i][0];
+    f[i][1] += qone * efield[i][1];
+    f[i][2] += qone * efield[i][2];
+    myeng += 0.5 * qone * epot[i];
+  } 
+
+  MPI_Allreduce(&myeng,&energy,1,MPI_DOUBLE,MPI_SUM,world);
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local arrays
+------------------------------------------------------------------------- */
+
+double Scafacos::memory_usage()
+{
+  double bytes = 0.0;
+  bytes += maxatom * sizeof(double);
+  bytes += 3*maxatom * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+    setup of ScaFaCoS handle with common parameters 
+------------------------------------------------------------------------- */
+
+void Scafacos::setup_handle()
+{
+  // store simulation box params
+  // NOTE: this assumes orthogonal box
+  // NOTE: total particles may not be a 4-byte int
+  // NOTE: what does SCFCS mean by offset?
+  //       it's an integer flag in LAMMPS
+
+  old_periodicity[0] = domain->xperiodic;
+  old_periodicity[1] = domain->yperiodic;
+  old_periodicity[2] = domain->zperiodic;
+
+  old_offset[0] = domain->boundary[0][0];
+  old_offset[1] = domain->boundary[1][0];
+  old_offset[2] = domain->boundary[2][0];
+
+  old_box_x[0] = domain->prd[0];
+  old_box_x[1] = old_box_x[2] = 0.0;
+  old_box_y[1] = domain->prd[1];
+  old_box_y[0] = old_box_y[2] = 0.0;
+  old_box_z[2] = domain->prd[2];
+  old_box_z[1] = old_box_z[0] = 0.0;
+
+  old_natoms = atom->natoms;
+
+  // set all required ScaFaCoS params
+
+  result = fcs_set_box_a(fcs,old_box_x);
+  check_result(result);
+
+  result = fcs_set_box_b(fcs,old_box_y);
+  check_result(result);
+
+  result = fcs_set_box_c(fcs,old_box_z);
+  check_result(result);
+
+  result = fcs_set_box_origin(fcs,old_offset);
+  check_result(result);
+
+  result = fcs_set_periodicity(fcs,old_periodicity);
+  check_result(result);
+
+  result = fcs_set_total_particles(fcs,old_natoms);
+  check_result(result);
+
+  // NOTE: for now disable short-range calculations within LAMMPS
+
+  int near_field_flag = 0;
+  result = fcs_set_near_field_flag(fcs,near_field_flag);
+  check_result(result);
+}
+
+/* ----------------------------------------------------------------------
+    check if box parameters changed, requiring a new call to fcs_tune
+------------------------------------------------------------------------- */
+
+bool Scafacos::box_has_changed()
+{
+  int *periodicity = domain->periodicity;
+  double *prd = domain->prd;
+
+  bool changed = 
+    (periodicity[0] != old_periodicity[0]) ||
+    (periodicity[1] != old_periodicity[1]) ||
+    (periodicity[2] != old_periodicity[2]) ||
+    (domain->boundary[0][0] != old_offset[0]) ||
+    (domain->boundary[1][0] != old_offset[1]) ||
+    (domain->boundary[2][0] != old_offset[2]) ||
+    (prd[0] != old_box_x[0]) ||
+    (prd[1] != old_box_y[1]) ||
+    (prd[2] != old_box_z[2]) ||
+    (atom->natoms != old_natoms);
+
+  return changed;
+}
+
+/* ----------------------------------------------------------------------
+   check ScaFaCoS result for error condition
+   NOTE: will all procs have same error?
+------------------------------------------------------------------------- */
+
+void Scafacos::check_result(FCSResult result) 
+{
+  if (!result) return;
+
+  std::stringstream ss;
+  ss << "ScaFaCoS: " << fcs_result_get_function(result) << "\n" 
+     << fcs_result_get_message(result) << "\n";
+  fcs_result_destroy(result);
+  std::string err_msg = ss.str();
+  const char *str = err_msg.c_str();
+
+  error->one(FLERR,str);
+}
+ 

src/USER-SCAFACOS/scafacos.h

+/* -*- 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(scafacos,Scafacos)
+
+#else
+
+#ifndef LMP_SCAFACOS_H
+#define LMP_SCAFACOS_H
+
+#include "kspace.h"
+#include "fcs.h"
+
+namespace LAMMPS_NS {
+
+class Scafacos : public KSpace {
+ public:
+  Scafacos(class LAMMPS *, int, char **);
+  ~Scafacos();
+  void init();
+  void setup();
+  void compute(int, int);
+  double memory_usage();
+
+ private: 
+  int me;
+
+  char *method;
+  double tolerance;
+  double *epot,**efield;
+  int tolerance_type;
+  int initialized,maxatom;
+
+  FCS fcs;                // ScaFaCoS handle
+  FCSResult result;       // result for each ScaFaCoS call
+
+  // simulation state: box, natoms
+  // so ScaFaCoS can detect if changes, e.g. for NPT
+
+  fcs_float old_box_x[3],old_box_y[3],old_box_z[3];
+  fcs_float old_offset[3];
+  fcs_int old_periodicity[3];
+  fcs_int old_natoms;
+
+  void check_result(FCSResult);
+  void setup_handle();
+  bool box_has_changed();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+*/