Permalink
Browse files

On-the-fly tuning of MMM1D

  • Loading branch information...
1 parent 9c78722 commit 79325e52a524c695fc5190a1ee9252d185387c56 @arnolda arnolda committed Feb 28, 2014
Showing with 97 additions and 112 deletions.
  1. +0 −1 src/interaction_data.cpp
  2. +90 −74 src/mmm1d.cpp
  3. +1 −16 src/mmm1d.hpp
  4. +2 −19 src/tcl/mmm1d_tcl.cpp
  5. +4 −2 testsuite/mmm1d.tcl
@@ -843,7 +843,6 @@ int coulomb_set_bjerrum(double bjerrum)
rf_params.B = 0.0;
case COULOMB_MMM1D:
mmm1d_params.maxPWerror = 1e40;
- mmm1d_params.bessel_cutoff = 0;
}
mpi_bcast_coulomb_params();
View
@@ -40,11 +40,12 @@
/** How many trial calculations */
#define TEST_INTEGRATIONS 1000
-/** Largest reasonable cutoff for Bessel function */
+/** Largest numerically stable cutoff for Bessel function. Don't
+ change without improving the formulas. */
#define MAXIMAL_B_CUT 30
-/** Granularity of the radius scan in multiples of box_l[2] */
-#define RAD_STEPPING 0.1
+/** minimal radius for the far formula in multiples of box_l[2] */
+#define MIN_RAD 0.01
/** if you define this, the Besselfunctions are calculated up
to machine precision, otherwise 10^-14, which should be
@@ -61,62 +62,65 @@
static double uz, L2, uz2, prefuz2, prefL3_i;
/*@}*/
-MMM1D_struct mmm1d_params = { 0.05, 5, 1, 1e-5 };
-
+MMM1D_struct mmm1d_params = { 0.05, 1e-5 };
+/** From which distance a certain bessel cutoff is valid. Can't be part of the
+ params since these get broadcasted. */
+static double *bessel_radii;
+
-static void MMM1D_setup_constants()
+static double far_error(int P, double minrad)
{
- uz = 1/box_l[2];
- L2 = box_l[2]*box_l[2];
- uz2 = uz*uz;
- prefuz2 = coulomb.prefactor*uz2;
- prefL3_i = prefuz2*uz;
+ // this uses an upper bound to all force components and the potential
+ double rhores = 2*M_PI*uz*minrad;
+ double pref = 4*uz*dmax(1, 2*M_PI*uz);
+
+ return pref*K1(rhores*P)*exp(rhores)/rhores*(P - 1 + 1/rhores);
}
-double determine_bessel_cutoff(double switch_rad, double maxPWerror, int maxP)
+static double determine_minrad(double maxPWerror, int P)
{
- /* this calculates an upper bound to all force components and the potential */
-
- // fprintf(stderr, "determ: %f %f %d\n", switch_rad, maxPWerror, maxP);
+ // bisection to search for where the error is maxPWerror
+ double rgranularity = MIN_RAD*box_l[2];
+ double rmin = rgranularity;
+ double rmax = dmin(box_l[0], box_l[1]);
+ double errmin = far_error(P, rmin);
+ double errmax = far_error(P, rmax);
+ if (errmin < maxPWerror) {
+ // we can do almost all radii with this P
+ return rmin;
+ }
+ if (errmax > maxPWerror) {
+ // make sure that this switching radius cannot be reached
+ return 2*dmax(box_l[0], box_l[1]);
+ }
- double err;
- double rhores = 2*M_PI*uz*switch_rad;
- double pref = 4*uz*dmax(1, 2*M_PI*uz);
- int P = 1;
- do {
- err = pref*K1(rhores*P)*exp(rhores)/rhores*(P - 1 + 1/rhores);
- // fprintf(stderr, "%d %e\n", P, err); */
- P++;
- } while (err > maxPWerror && P <= maxP);
- P--;
- return P;
+ while (rmax - rmin > rgranularity) {
+ double c = 0.5*(rmin + rmax);
+ double errc = far_error(P, c);
+ if (errc > maxPWerror) {
+ rmin = c;
+ } else {
+ rmax = c;
+ }
+ }
+ return 0.5*(rmin + rmax);
}
-int MMM1D_set_params(double switch_rad, int bessel_cutoff, double maxPWerror)
+static void determine_bessel_radii(double maxPWerror, int maxP)
{
- MMM1D_setup_constants();
-
- mmm1d_params.far_switch_radius_2 = (switch_rad > 0) ? SQR(switch_rad) : -1;
- mmm1d_params.bessel_cutoff = bessel_cutoff;
- /* if parameters come from here they are never calculated that is
- only the case if you call mmm1d_tune, which then
- changes this flag */
- mmm1d_params.bessel_calculated = 0;
-
- mmm1d_params.maxPWerror = maxPWerror;
- coulomb.method = COULOMB_MMM1D;
-
- mpi_bcast_coulomb_params();
-
- return 0;
+ bessel_radii = (double *)realloc(bessel_radii, sizeof(double)*maxP);
+ for (int P = 1; P <= maxP; ++P) {
+ bessel_radii[P-1] = determine_minrad(maxPWerror, P);
+ //printf("cutoff %d %f\n", P, bessel_radii[P-1]);
+ }
}
-void MMM1D_recalcTables()
+static void prepare_polygamma_series(double maxPWerror, double maxrad2)
{
/* polygamma, determine order */
int n;
double err;
- double rhomax2nm2, rhomax2 = uz2*mmm1d_params.far_switch_radius_2;
+ double rhomax2nm2, rhomax2 = uz2*maxrad2;
/* rhomax2 < 1, so rhomax2m2 falls monotonously */
n = 1;
@@ -130,7 +134,18 @@ void MMM1D_recalcTables()
n++;
// fprintf(stderr, "%f\n", err);
}
- while (err > 0.1*mmm1d_params.maxPWerror);
+ while (err > 0.1*maxPWerror);
+}
+
+int MMM1D_set_params(double switch_rad, double maxPWerror)
+{
+ mmm1d_params.far_switch_radius_2 = (switch_rad > 0) ? SQR(switch_rad) : -1;
+ mmm1d_params.maxPWerror = maxPWerror;
+ coulomb.method = COULOMB_MMM1D;
+
+ mpi_bcast_coulomb_params();
+
+ return 0;
}
int MMM1D_sanity_checks()
@@ -157,12 +172,14 @@ void MMM1D_init()
if (mmm1d_params.far_switch_radius_2 >= SQR(box_l[2]))
mmm1d_params.far_switch_radius_2 = 0.8*SQR(box_l[2]);
- MMM1D_setup_constants();
+ uz = 1/box_l[2];
+ L2 = box_l[2]*box_l[2];
+ uz2 = uz*uz;
+ prefuz2 = coulomb.prefactor*uz2;
+ prefL3_i = prefuz2*uz;
- if (mmm1d_params.bessel_calculated) {
- mmm1d_params.bessel_cutoff = determine_bessel_cutoff(sqrt(mmm1d_params.far_switch_radius_2), mmm1d_params.maxPWerror, MAXIMAL_B_CUT);
- }
- MMM1D_recalcTables();
+ determine_bessel_radii(mmm1d_params.maxPWerror, MAXIMAL_B_CUT);
+ prepare_polygamma_series(mmm1d_params.maxPWerror, mmm1d_params.far_switch_radius_2);
}
void add_mmm1d_coulomb_pair_force(double chpref, double d[3], double r2, double r, double force[3])
@@ -242,7 +259,10 @@ void add_mmm1d_coulomb_pair_force(double chpref, double d[3], double r2, double
double sr = 0, sz = 0;
int bp;
- for (bp = 1; bp < mmm1d_params.bessel_cutoff; bp++) {
+ for (bp = 1; bp < MAXIMAL_B_CUT; bp++) {
+ if (bessel_radii[bp-1] < rxy)
+ break;
+
double fq = C_2PI*bp, k0, k1;
#ifdef BESSEL_MACHINE_PREC
k0 = K0(fq*rxy_d);
@@ -320,7 +340,10 @@ double mmm1d_coulomb_pair_energy(Particle *p1, Particle *p2, double d[3], double
/* The first Bessel term will compensate a little bit the
log term, so add them close together */
E = -0.25*log(rxy2_d) + 0.5*(M_LN2 - C_GAMMA);
- for (bp = 1; bp < mmm1d_params.bessel_cutoff; bp++) {
+ for (bp = 1; bp < MAXIMAL_B_CUT; bp++) {
+ if (bessel_radii[bp-1] < rxy)
+ break;
+
double fq = C_2PI*bp;
E += K0(fq*rxy_d)*cos(fq*z_d);
}
@@ -337,15 +360,16 @@ int mmm1d_tune(char **log)
double maxrad = box_l[2]; /* N_psi = 2, theta=2/3 maximum for rho */
double switch_radius;
- if (mmm1d_params.bessel_cutoff < 0 && mmm1d_params.far_switch_radius_2 < 0) {
- /* determine besselcutoff and optimal switching radius */
- for (switch_radius = RAD_STEPPING*maxrad;
- switch_radius < maxrad;
- switch_radius += RAD_STEPPING*maxrad) {
- mmm1d_params.bessel_cutoff = determine_bessel_cutoff(switch_radius, mmm1d_params.maxPWerror, MAXIMAL_B_CUT);
- /* no reasonable cutoff possible */
- if (mmm1d_params.bessel_cutoff == MAXIMAL_B_CUT)
- continue;
+ if (mmm1d_params.far_switch_radius_2 < 0) {
+ /* determine besselcutoff and optimal switching radius. Should be around 0.33 */
+ for (switch_radius = 0.2*maxrad;
+ switch_radius < 0.4*maxrad;
+ switch_radius += 0.025*maxrad) {
+ if (switch_radius <= bessel_radii[MAXIMAL_B_CUT - 1]) {
+ // this switching radius is too small for our Bessel series
+ continue;
+ }
+
mmm1d_params.far_switch_radius_2 = SQR(switch_radius);
coulomb.method = COULOMB_MMM1D;
@@ -360,8 +384,7 @@ int mmm1d_tune(char **log)
if (int_time < 0)
return ES_ERROR;
- sprintf(buffer, "r= %f c= %d t= %f ms\n",
- switch_radius, mmm1d_params.bessel_cutoff, int_time);
+ sprintf(buffer, "r= %f t= %f ms\n", switch_radius, int_time);
*log = strcat_alloc(*log, buffer);
if (int_time < min_time) {
@@ -372,23 +395,16 @@ int mmm1d_tune(char **log)
else if (int_time > 2*min_time)
break;
}
- switch_radius = min_rad;
+ switch_radius = min_rad;
mmm1d_params.far_switch_radius_2 = SQR(switch_radius);
- mmm1d_params.bessel_cutoff = determine_bessel_cutoff(switch_radius, mmm1d_params.maxPWerror, MAXIMAL_B_CUT);
- mmm1d_params.bessel_calculated = 1;
}
- else if (mmm1d_params.bessel_cutoff < 0) {
- /* determine besselcutoff to achieve at least the given pairwise error */
- mmm1d_params.bessel_cutoff = determine_bessel_cutoff(sqrt(mmm1d_params.far_switch_radius_2),
- mmm1d_params.maxPWerror, MAXIMAL_B_CUT);
- if (mmm1d_params.bessel_cutoff == MAXIMAL_B_CUT) {
+ else {
+ if (mmm1d_params.far_switch_radius_2 <= SQR(bessel_radii[MAXIMAL_B_CUT - 1])) {
+ // this switching radius is too small for our Bessel series
*log = strcat_alloc(*log, "could not find reasonable bessel cutoff");
return ES_ERROR;
}
- mmm1d_params.bessel_calculated = 1;
}
- else
- mmm1d_params.bessel_calculated = 0;
coulomb.method = COULOMB_MMM1D;
View
@@ -38,16 +38,6 @@
typedef struct {
/** square of the switching radius */
double far_switch_radius_2;
- /** cutoff of the bessel sum */
- int bessel_cutoff;
- /** Wether to recalculate the Bessel cutoff automatically.
- If some parameters like the box dimensions change, the current
- Bessel cutoff may not be suitable to achieve the required accuracy
- anymore. If the user did not specify the Bessel cutoff explicitly,
- this flag is 1, and whenever necessary, the Bessel cutoff is
- recalculated.
- */
- int bessel_calculated;
/** required accuracy */
double maxPWerror;
} MMM1D_struct;
@@ -58,14 +48,9 @@ extern MMM1D_struct mmm1d_params;
if you set this parameters.
@param switch_rad at which xy-distance the calculation switches from the far to the
near formula. If -1, this parameter will be tuned automatically.
- @param bessel_cutoff the cutoff for the bessel sum, aka far formula. Normally set this
- to -1, then the cutoff is automatically determined using the error formula.
@param maxPWerror the maximal allowed error for the potential and the forces without the
prefactors, i. e. for the pure lattice 1/r-sum. */
-int MMM1D_set_params(double switch_rad, int bessel_cutoff, double maxPWerror);
-
-/** recalculate the polygamma taylor series. */
-void MMM1D_recalcTables();
+int MMM1D_set_params(double switch_rad, double maxPWerror);
/// check that MMM1D can run with the current parameters
int MMM1D_sanity_checks();
View
@@ -46,8 +46,6 @@ int tclprint_to_result_MMM1D(Tcl_Interp *interp)
Tcl_PrintDouble(interp, sqrt(mmm1d_params.far_switch_radius_2), buffer);
Tcl_AppendResult(interp, "mmm1d ", buffer, " ",(char *) NULL);
- sprintf(buffer, "%d", mmm1d_params.bessel_cutoff);
- Tcl_AppendResult(interp, buffer, " ",(char *) NULL);
Tcl_PrintDouble(interp, mmm1d_params.maxPWerror, buffer);
Tcl_AppendResult(interp, buffer,(char *) NULL);
@@ -57,7 +55,6 @@ int tclprint_to_result_MMM1D(Tcl_Interp *interp)
int tclcommand_inter_coulomb_parse_mmm1d(Tcl_Interp *interp, int argc, char **argv)
{
double switch_rad, maxPWerror;
- int bessel_cutoff;
if (argc < 2) {
Tcl_AppendResult(interp, "wrong # arguments: inter coulomb mmm1d <switch radius> "
@@ -69,7 +66,6 @@ int tclcommand_inter_coulomb_parse_mmm1d(Tcl_Interp *interp, int argc, char **ar
/* autodetermine bessel cutoff AND switching radius */
if (! ARG_IS_D(1, maxPWerror))
return TCL_ERROR;
- bessel_cutoff = -1;
switch_rad = -1;
}
else {
@@ -78,23 +74,10 @@ int tclcommand_inter_coulomb_parse_mmm1d(Tcl_Interp *interp, int argc, char **ar
if ((! ARG_IS_D(0, switch_rad)) ||
(! ARG_IS_D(1, maxPWerror)))
return TCL_ERROR;
- bessel_cutoff = -1;
- }
- else if (argc == 3) {
- /* fully manual */
- if((! ARG_IS_D(0, switch_rad)) ||
- (! ARG_IS_I(1, bessel_cutoff)) ||
- (! ARG_IS_D(2, maxPWerror)))
- return TCL_ERROR;
-
- if (bessel_cutoff <=0) {
- Tcl_AppendResult(interp, "bessel cutoff too small", (char *)NULL);
- return TCL_ERROR;
- }
}
else {
Tcl_AppendResult(interp, "wrong # arguments: inter coulomb mmm1d <switch radius> "
- "{<bessel cutoff>} <maximal error for near formula> | tune <maximal pairwise error>", (char *) NULL);
+ "<maximal error for near formula> | tune <maximal pairwise error>", (char *) NULL);
return TCL_ERROR;
}
@@ -104,7 +87,7 @@ int tclcommand_inter_coulomb_parse_mmm1d(Tcl_Interp *interp, int argc, char **ar
}
}
- MMM1D_set_params(switch_rad, bessel_cutoff, maxPWerror);
+ MMM1D_set_params(switch_rad, maxPWerror);
char *log = NULL;
int result = mmm1d_tune(&log) == ES_OK ? TCL_OK : TCL_ERROR;
View
@@ -57,7 +57,7 @@ if { [catch {
############## mmm1d-specific part
setmd periodic 0 0 1
- inter coulomb 1.0 mmm1d 6.0 3 0.0001
+ inter coulomb 1.0 mmm1d 6.0 0.0001
# to ensure force recalculation
invalidate_system
@@ -66,7 +66,9 @@ if { [catch {
# here you can create the necessary snapshot
if { 0 } {
inter coulomb 1.0 mmm1d tune 1e-20
- write_data "mmm1d_system.data22"
+ puts [inter coulomb]
+ integrate 0
+ write_data "mmm1d_system.data"
}
############## end

0 comments on commit 79325e5

Please sign in to comment.