diff --git a/src/USER-SCAFACOS/scafacos.cpp b/src/USER-SCAFACOS/scafacos.cpp
index fef4b4948fde404a8b4a1eb8acf3f391dde36393..8779796a5431dc9fc41df52492e8bfa172bb4f75 100644
--- a/src/USER-SCAFACOS/scafacos.cpp
+++ b/src/USER-SCAFACOS/scafacos.cpp
@@ -38,16 +38,24 @@ 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");
+ if (narg > 2) error->all(FLERR,"Illegal scafacos command");
int n = strlen(arg[0]) + 1;
method = new char[n];
strcpy(method,arg[0]);
- tolerance = force->numeric(FLERR,arg[1]);
+ // to allow 'kspace scafacos <method>' with default values
+ if (narg == 2)
+ tolerance = force->numeric(FLERR,arg[1]);
+ else
+ tolerance = 0.001;
// optional ScaFaCoS library setting defaults
-
- tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
+ // choose the correct default tolerance type for chosen method
+ // TODO: needs to be expanded for all solvers, currently mainly used ones
+ if (strcmp(method,"fmm") == 0)
+ tolerance_type = FCS_TOLERANCE_TYPE_ENERGY;
+ else if (strcmp(method,"p3m") == 0 || strcmp(method,"p2nfft") == 0)
+ tolerance_type = FCS_TOLERANCE_TYPE_FIELD;
// initializations
@@ -68,7 +76,8 @@ Scafacos::~Scafacos()
memory->destroy(epot);
memory->destroy(efield);
- // RENE: any clean-up/shut-down call to ScaFaCoS needed?
+ // clean up of the ScaFaCoS handle and internal arrays
+ fcs_destroy(fcs);
}
/* ---------------------------------------------------------------------- */
@@ -93,7 +102,6 @@ void Scafacos::init()
// one-time initialization of ScaFaCoS
- scale = 1.0;
qqrd2e = force->qqrd2e;
if (!initialized) {
@@ -129,9 +137,7 @@ void Scafacos::compute(int eflag, int vflag)
double *q = atom->q;
int nlocal = atom->nlocal;
- // RENE: why is scale needed?
-
- const double qscale = qqrd2e * scale;
+ const double qscale = qqrd2e;
if (eflag || vflag) ev_setup(eflag,vflag);
else eflag_atom = 0;
@@ -155,16 +161,6 @@ void Scafacos::compute(int eflag, int vflag)
}
- // initialize epot & efield
- // RENE: is this necessary? or does Scafacos just set them
-
- 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);
@@ -196,7 +192,7 @@ void Scafacos::compute(int eflag, int vflag)
int Scafacos::modify_param(int narg, char **arg)
{
- // RENE: add any Scafacos options here you want to expose to LAMMPS
+ // add any Scafacos options here you want to expose to LAMMPS
// syntax: kspace_modify scafacos keyword value1 value2 ...
// keyword = tolerance
// value1 = energy, energy_rel, etc
@@ -245,17 +241,17 @@ void Scafacos::setup_handle()
{
// store simulation box params
+ // setup periodicity
old_periodicity[0] = domain->xperiodic;
old_periodicity[1] = domain->yperiodic;
old_periodicity[2] = domain->zperiodic;
- // RENE: what does SCFCS mean by offset?
- // it's an integer flag in LAMMPS, but being stored in a float?
-
- old_offset[0] = domain->boundary[0][0];
- old_offset[1] = domain->boundary[1][0];
- old_offset[2] = domain->boundary[2][0];
+ // setup box origin (lower left front corner of the system)
+ old_origin[0] = domain->boundary[0][0];
+ old_origin[1] = domain->boundary[1][0];
+ old_origin[2] = domain->boundary[2][0];
+ // setup box vectors (base vectors of the system box)
old_box_x[0] = domain->prd[0];
old_box_x[1] = old_box_x[2] = 0.0;
old_box_y[1] = domain->prd[1];
@@ -263,10 +259,10 @@ void Scafacos::setup_handle()
old_box_z[2] = domain->prd[2];
old_box_z[1] = old_box_z[0] = 0.0;
+ // setup number of atoms in the system
old_natoms = atom->natoms;
- // set all required ScaFaCoS params
-
+ // store parameters to ScaFaCoS handle
result = fcs_set_box_a(fcs,old_box_x);
check_result(result);
@@ -276,7 +272,7 @@ void Scafacos::setup_handle()
result = fcs_set_box_c(fcs,old_box_z);
check_result(result);
- result = fcs_set_box_origin(fcs,old_offset);
+ result = fcs_set_box_origin(fcs,old_origin);
check_result(result);
result = fcs_set_periodicity(fcs,old_periodicity);
@@ -285,11 +281,12 @@ void Scafacos::setup_handle()
result = fcs_set_total_particles(fcs,old_natoms);
check_result(result);
- // RENE: disable short-range calculations within LAMMPS
- // not sure what this is doing
- // is this the correct thing to do for now?
-
- int near_field_flag = 0;
+ // allow ScaFaCoS to calculate the near field computations for now
+ // TODO: allow the delegation of the near field computations
+ // within LAMMPS
+ // (near_field_flag = 1 -> enables the internal near field calcs
+ // 0 -> disables the internal near field calcs
+ int near_field_flag = 1;
result = fcs_set_near_field_flag(fcs,near_field_flag);
check_result(result);
}
@@ -307,9 +304,9 @@ bool Scafacos::box_has_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]) ||
+ (domain->boundary[0][0] != old_origin[0]) ||
+ (domain->boundary[1][0] != old_origin[1]) ||
+ (domain->boundary[2][0] != old_origin[2]) ||
(prd[0] != old_box_x[0]) ||
(prd[1] != old_box_y[1]) ||
(prd[2] != old_box_z[2]) ||
@@ -333,9 +330,6 @@ void Scafacos::check_result(FCSResult result)
std::string err_msg = ss.str();
const char *str = err_msg.c_str();
- // RENE: will all procs have same error?
- // if so, then should call error->all(FLERR,str)
-
error->one(FLERR,str);
}
diff --git a/src/USER-SCAFACOS/scafacos.h b/src/USER-SCAFACOS/scafacos.h
index f3708129f1c619f7a1926891bc647473a31922a0..3e549cf64ee1ba512f3b697161091d1fb62249bb 100644
--- a/src/USER-SCAFACOS/scafacos.h
+++ b/src/USER-SCAFACOS/scafacos.h
@@ -51,7 +51,7 @@ class Scafacos : public KSpace {
// 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_float old_origin[3];
fcs_int old_periodicity[3];
fcs_int old_natoms;