msmpot_setup.c

00001 /***************************************************************************
00002 *cr
00003 *cr (C) Copyright 2008-2009 The Board of Trustees of the
00004 *cr University of Illinois
00005 *cr All Rights Reserved
00006 *cr
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 * $RCSfile: msmpot_setup.c,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.7 $ $Date: 2020年10月21日 20:35:05 $
00015 *
00016 ***************************************************************************/
00017 
00018 #include "msmpot_internal.h"
00019 
00020 
00021 int Msmpot_check_params(Msmpot *msm, const float *epotmap,
00022 int mx, int my, int mz, float lx, float ly, float lz,
00023 float vx, float vy, float vz, const float *atom, int natoms) {
00024 
00025 if (NULL == epotmap) {
00026 return ERRMSG(MSMPOT_ERROR_PARAM, "map buffer is NULL");
00027 }
00028 else if (NULL == atom) {
00029 return ERRMSG(MSMPOT_ERROR_PARAM, "atom array is NULL");
00030 }
00031 else if (natoms <= 0) {
00032 return ERRMSG(MSMPOT_ERROR_PARAM, "number of atoms is not positive");
00033 }
00034 else if (mx <= 0 || my <= 0 || mz <= 0) {
00035 return ERRMSG(MSMPOT_ERROR_PARAM, "number of map points is not positive");
00036 }
00037 else if (lx <= 0 || ly <= 0 || lz <= 0) {
00038 return ERRMSG(MSMPOT_ERROR_PARAM, "map lengths are not positive");
00039 }
00040 else if ((vx <= msm->hmin && vx != 0) || (vy <= msm->hmin && vy != 0)
00041 || (vz <= msm->hmin && vz != 0)) {
00042 return ERRMSG(MSMPOT_ERROR_PARAM,
00043 "periodic cell lengths must be greater than MSM spacing");
00044 }
00045 else if ((lx > vx && vx != 0) || (ly > vy && vy != 0)
00046 || (lz > vz && vz != 0)) {
00047 return ERRMSG(MSMPOT_ERROR_PARAM,
00048 "map lengths must be within periodic cell");
00049 }
00050 
00051 return MSMPOT_SUCCESS;
00052 }
00053 
00054 
00055 /* called by Msmpot_create() */
00056 void Msmpot_set_defaults(Msmpot *msm) {
00057 /* setup default parameters */
00058 msm->hmin = DEFAULT_HMIN;
00059 msm->a = DEFAULT_CUTOFF;
00060 msm->interp = DEFAULT_INTERP;
00061 msm->split = DEFAULT_SPLIT;
00062 msm->bindepth = DEFAULT_BINDEPTH;
00063 msm->binfill = DEFAULT_BINFILL;
00064 msm->density = DEFAULT_DENSITY;
00065 msm->errtol = ((float) DEFAULT_ERRTOL);
00066 
00067 #ifdef MSMPOT_CUDA
00068 /* attempt to use CUDA but fall back on CPU if unable */
00069 msm->use_cuda = 1;
00070 msm->cuda_optional = 1;
00071 #endif
00072 }
00073 
00074 
00075 int Msmpot_configure(Msmpot *msm,
00076 int interp,
00077 int split,
00078 float cutoff,
00079 float hmin,
00080 int nlevels,
00081 float density,
00082 float binfill,
00083 float errtol,
00084 int usecuda
00085 ) {
00086 if (interp < 0 || interp >= MSMPOT_INTERPMAX) {
00087 return ERROR(MSMPOT_ERROR_PARAM);
00088 }
00089 else msm->interp = interp; /* default is 0 */
00090 
00091 if (split < 0 || split >= MSMPOT_SPLITMAX) {
00092 return ERROR(MSMPOT_ERROR_PARAM);
00093 }
00094 else msm->split = split; /* default is 0 */
00095 
00096 if (cutoff < 0) return ERROR(MSMPOT_ERROR_PARAM);
00097 else if (cutoff > 0) msm->a = cutoff;
00098 else msm->a = DEFAULT_CUTOFF;
00099 
00100 if (hmin < 0) return ERROR(MSMPOT_ERROR_PARAM);
00101 else if (hmin > 0) msm->hmin = hmin;
00102 else msm->hmin = DEFAULT_HMIN;
00103 
00104 if (nlevels < 0) return ERROR(MSMPOT_ERROR_PARAM);
00105 else msm->nlevels = nlevels; /* 0 is default */
00106 
00107 if (density < 0) return ERROR(MSMPOT_ERROR_PARAM);
00108 else if (density > 0) msm->density = density;
00109 else msm->density = DEFAULT_DENSITY;
00110 
00111 if (binfill < 0) return ERROR(MSMPOT_ERROR_PARAM);
00112 else if (binfill > 0) msm->binfill = binfill;
00113 else msm->binfill = DEFAULT_BINFILL;
00114 
00115 if (errtol < 0) return ERROR(MSMPOT_ERROR_PARAM);
00116 else if (errtol > 0) msm->errtol = errtol;
00117 else msm->errtol = ((float) DEFAULT_ERRTOL);
00118 
00119 #ifdef MSMPOT_CUDA
00120 msm->use_cuda = (usecuda != 0);
00121 #else
00122 if (usecuda != 0) {
00123 return ERROR(MSMPOT_ERROR_SUPPORT);
00124 }
00125 #endif
00126 
00127 return MSMPOT_SUCCESS;
00128 }
00129 
00130 
00131 /* called by Msmpot_destroy() */
00132 void Msmpot_cleanup(Msmpot *msm) {
00133 int i;
00134 for (i = 0; i < msm->maxlevels; i++) {
00135 GRID_DONE( &(msm->qh[i]) );
00136 GRID_DONE( &(msm->eh[i]) );
00137 GRID_DONE( &(msm->gc[i]) );
00138 }
00139 free(msm->qh);
00140 free(msm->eh);
00141 free(msm->gc);
00142 free(msm->ezd);
00143 free(msm->eyzd);
00144 free(msm->lzd);
00145 free(msm->lyzd);
00146 free(msm->phi_x);
00147 free(msm->phi_y);
00148 free(msm->phi_z);
00149 
00150 free(msm->first_atom_index);
00151 free(msm->next_atom_index);
00152 
00153 free(msm->bin);
00154 free(msm->bincount);
00155 free(msm->over);
00156 free(msm->boff);
00157 }
00158 
00159 
00160 static int setup_domain(Msmpot *msm);
00161 static int setup_bins(Msmpot *msm);
00162 static int setup_origin(Msmpot *msm);
00163 
00164 static int setup_hierarchy(Msmpot *msm);
00165 
00166 /* called by setup_hierarchy() */
00167 static int setup_periodic_hlevelparams_1d(Msmpot *msm,
00168 float len, /* domain length */
00169 float *hh, /* determine h */
00170 int *nn, /* determine number grid spacings covering domain */
00171 int *aindex, /* determine smallest lattice index */
00172 int *bindex /* determine largest lattice index */
00173 );
00174 
00175 /* called by setup_hierarchy() */
00176 static int setup_nonperiodic_hlevelparams_1d(Msmpot *msm,
00177 float len, /* measure to furthest point interpolated from grid */
00178 float *hh, /* determine h */
00179 int *nn, /* determine number grid spacings covering domain */
00180 int *aindex, /* determine smallest lattice index */
00181 int *bindex /* determine largest lattice index */
00182 );
00183 
00184 static int setup_mapinterp(Msmpot *msm);
00185 
00186 static int setup_mapinterpcoef_1d(Msmpot *msm,
00187 float h, /* spacing of MSM h-level lattice */
00188 float delta, /* spacing of epotmap lattice */
00189 int n, /* number of MSM h spacings to cover domain */
00190 int m, /* number of epotmap lattice points */
00191 float *p_h_delta, /* calculate ratio h/delta */
00192 int *p_cycle, /* number of MSM points until next map alignment */
00193 int *p_rmap, /* radius of map points about MSM point */
00194 float **p_phi, /* coefficients that weight map about MSM points */
00195 int *p_max_phi /* size of phi memory allocation */
00196 );
00197 
00198 
00199 #ifdef MSMPOT_VERBOSE
00200 static int print_status(Msmpot *msm);
00201 #endif
00202 
00203 
00204 /* called by Msmpot_compute() */
00205 int Msmpot_setup(Msmpot *msm) {
00206 int err = 0;
00207 
00208 REPORT("Setting up for MSM computation.");
00209 
00210 /* set domain lengths, needed for any non-periodic dimensions */
00211 err = setup_domain(msm);
00212 if (err) return ERROR(err);
00213 
00214 /* determine bin subdivision across domain */
00215 err = setup_bins(msm);
00216 if (err) return ERROR(err);
00217 
00218 /* given bin subdivision, find good domain origin for periodic dimensions */
00219 err = setup_origin(msm);
00220 if (err) return ERROR(err);
00221 
00222 /* set up hierarchy of lattices for long-range part */
00223 err = setup_hierarchy(msm);
00224 if (err) return ERROR(err);
00225 
00226 
00227 #if ! defined(MSMPOT_SHORTRANGE_ONLY)
00228 if (msm->px == msm->lx && msm->py == msm->ly && msm->pz == msm->lz) {
00229 /* determine map interpolation parameters
00230 * and MSM lattice spacings hx, hy, hz */
00231 err = setup_mapinterp(msm);
00232 if (err) return ERROR(err);
00233 }
00234 #endif
00235 
00236 
00237 #ifdef MSMPOT_VERBOSE
00238 err = print_status(msm);
00239 if (err) return ERROR(err);
00240 #endif
00241 
00242 #ifdef MSMPOT_CUDA
00243 /* set up CUDA device */
00244 if (msm->use_cuda) {
00245 err = Msmpot_cuda_setup(msm->msmcuda, msm);
00246 if (err) return ERROR(err);
00247 }
00248 #endif
00249 
00250 return MSMPOT_SUCCESS;
00251 }
00252 
00253 
00254 typedef struct InterpParams_t {
00255 int nu;
00256 int stencil;
00257 int omega;
00258 } InterpParams;
00259 
00260 static InterpParams INTERP_PARAMS[] = {
00261 { 1, 4, 6 }, /* cubic */
00262 { 2, 6, 10 }, /* quintic */
00263 { 2, 6, 10 }, /* quintic, C2 */
00264 { 3, 8, 14 }, /* septic */
00265 { 3, 8, 14 }, /* septic, C3 */
00266 { 4, 10, 18 }, /* nonic */
00267 { 4, 10, 18 }, /* nonic, C4 */
00268 };
00269 
00270 
00271 #ifdef MSMPOT_VERBOSE
00272 int print_status(Msmpot *msm) {
00273 int j, k;
00274 int ispx = (IS_SET_X(msm->isperiodic) != 0);
00275 int ispy = (IS_SET_Y(msm->isperiodic) != 0);
00276 int ispz = (IS_SET_Z(msm->isperiodic) != 0);
00277 int ispany = (IS_SET_ANY(msm->isperiodic) != 0);
00278 printf("#MSMPOT STATUS\n");
00279 printf("# natoms= %d\n", msm->natoms);
00280 printf("# domain lengths= %g %g %g\n", msm->px, msm->py, msm->pz);
00281 printf("# domain origin= %g %g %g\n", msm->px0, msm->py0, msm->pz0);
00282 printf("# map size= %d x %d x %d\n", msm->mx, msm->my, msm->mz);
00283 printf("# map spacing= %g, %g, %g\n", msm->dx, msm->dy, msm->dz);
00284 printf("# map lengths= %g, %g, %g\n", msm->lx, msm->ly, msm->lz);
00285 printf("# map origin= %g, %g, %g\n", msm->lx0, msm->ly0, msm->lz0);
00286 printf("# hmin= %g\n", msm->hmin);
00287 printf("# cutoff= %g\n", msm->a);
00288 printf("# MSM size= %d x %d x %d\n", msm->nx, msm->ny, msm->nz);
00289 printf("# MSM spacing= %g, %g, %g\n", msm->hx, msm->hy, msm->hz);
00290 printf("# MSM interpolation= %d\n", msm->interp);
00291 printf("# MSM splitting= %d\n", msm->split);
00292 printf("# MSM number of levels= %d\n", msm->nlevels);
00293 printf("# MSM lattice hierarchy:\n");
00294 for (k = 0; k < msm->nlevels; k++) {
00295 floatGrid *qh = &(msm->qh[k]);
00296 int ia = qh->i0;
00297 int ib = ia + qh->ni - 1;
00298 int ja = qh->j0;
00299 int jb = ja + qh->nj - 1;
00300 int ka = qh->k0;
00301 int kb = ka + qh->nk - 1;
00302 printf("# level= %d: [%d..%d] x [%d..%d] x [%d..%d]\n",
00303 k, ia, ib, ja, jb, ka, kb);
00304 }
00305 printf("# ispx= %d ispy= %d ispz= %d ispany= %d\n",
00306 ispx, ispy, ispz, ispany);
00307 printf("# hx_dx= %g hy_dy= %g hz_dz= %g\n",
00308 msm->hx_dx, msm->hy_dy, msm->hz_dz);
00309 printf("# cycle_x= %d rmap_x= %d\n", msm->cycle_x, msm->rmap_x);
00310 for (k = 0; k < msm->cycle_x; k++) {
00311 int jtotal = 2*msm->rmap_x + 1;
00312 float *phi = msm->phi_x + k*jtotal;
00313 float phisum = 0;
00314 for (j = 0; j < jtotal; j++) phisum += phi[j];
00315 printf("# %d: sum= %g (", k, phisum);
00316 for (j = 0; j < jtotal; j++) {
00317 printf("%s%g", (0==j ? "= " : " + "), phi[j]);
00318 }
00319 printf(")\n");
00320 }
00321 printf("# cycle_y= %d rmap_y= %d\n", msm->cycle_y, msm->rmap_y);
00322 for (k = 0; k < msm->cycle_y; k++) {
00323 int jtotal = 2*msm->rmap_y + 1;
00324 float *phi = msm->phi_y + k*jtotal;
00325 float phisum = 0;
00326 for (j = 0; j < jtotal; j++) phisum += phi[j];
00327 printf("# %d: sum= %g (", k, phisum);
00328 for (j = 0; j < jtotal; j++) {
00329 printf("%s%g", (0==j ? "= " : " + "), phi[j]);
00330 }
00331 printf(")\n");
00332 }
00333 printf("# cycle_z= %d rmap_z= %d\n", msm->cycle_z, msm->rmap_z);
00334 for (k = 0; k < msm->cycle_z; k++) {
00335 int jtotal = 2*msm->rmap_z + 1;
00336 float *phi = msm->phi_z + k*jtotal;
00337 float phisum = 0;
00338 for (j = 0; j < jtotal; j++) phisum += phi[j];
00339 printf("# %d: sum= %g (", k, phisum);
00340 for (j = 0; j < jtotal; j++) {
00341 printf("%s%g", (0==j ? "= " : " + "), phi[j]);
00342 }
00343 printf(")\n");
00344 }
00345 printf("# bin size= %g %g %g\n", msm->bx, msm->by, msm->bz);
00346 printf("# atom bins= %d x %d x %d\n",
00347 msm->nbx, msm->nby, msm->nbz);
00348 return MSMPOT_SUCCESS;
00349 }
00350 #endif
00351 
00352 
00353 int setup_mapinterp(Msmpot *msm) {
00354 int mymz = msm->my * msm->mz;
00355 int err = 0;
00356 
00357 ASSERT(msm->mx > 0);
00358 ASSERT(msm->my > 0);
00359 ASSERT(msm->mz > 0);
00360 if (msm->max_eyzd < mymz) {
00361 float *t;
00362 t = (float *) realloc(msm->eyzd, mymz * sizeof(float));
00363 if (NULL == t) return ERROR(MSMPOT_ERROR_MALLOC);
00364 msm->eyzd = t;
00365 msm->max_eyzd = mymz;
00366 }
00367 if (msm->max_ezd < msm->mz) {
00368 float *t;
00369 t = (float *) realloc(msm->ezd, msm->mz * sizeof(float));
00370 if (NULL == t) return ERROR(MSMPOT_ERROR_MALLOC);
00371 msm->ezd = t;
00372 msm->max_ezd = msm->mz;
00373 }
00374 
00375 if (msm->px != msm->lx
00376 || msm->py != msm->ly
00377 || msm->pz != msm->lz) {
00378 /* must bail, can't do interpolation right yet */
00379 printf("px = %f lx = %f\n", msm->px, msm->lx);
00380 return ERRMSG(MSMPOT_ERROR_SUPPORT,
00381 "can't do interpolation for map lengths "
00382 "not equal to atom domain lengths");
00383 }
00384 
00385 
00386 err |= setup_mapinterpcoef_1d(msm,
00387 msm->hx, msm->dx, msm->nx, msm->mx, &(msm->hx_dx),
00388 &(msm->cycle_x), &(msm->rmap_x), &(msm->phi_x), &(msm->max_phi_x));
00389 err |= setup_mapinterpcoef_1d(msm,
00390 msm->hy, msm->dy, msm->ny, msm->my, &(msm->hy_dy),
00391 &(msm->cycle_y), &(msm->rmap_y), &(msm->phi_y), &(msm->max_phi_y));
00392 err |= setup_mapinterpcoef_1d(msm,
00393 msm->hz, msm->dz, msm->nz, msm->mz, &(msm->hz_dz),
00394 &(msm->cycle_z), &(msm->rmap_z), &(msm->phi_z), &(msm->max_phi_z));
00395 if (err) return ERROR(err);
00396 
00397 
00398 return MSMPOT_SUCCESS;
00399 }
00400 
00401 
00402 static int gcd(int a, int b) {
00403 /* subtraction-based Euclidean algorithm, from Knuth */
00404 if (0 == a) return b;
00405 while (b != 0) {
00406 if (a > b) a -= b;
00407 else b -= a;
00408 }
00409 return a;
00410 }
00411 
00412 
00413 int setup_mapinterpcoef_1d(Msmpot *msm,
00414 float h, /* spacing of MSM h-level lattice */
00415 float delta, /* spacing of epotmap lattice */
00416 int n, /* number of MSM h spacings to cover domain */
00417 int m, /* number of epotmap lattice points */
00418 float *p_h_delta, /* calculate ratio h/delta */
00419 int *p_cycle, /* number of MSM points until next map alignment */
00420 int *p_rmap, /* radius of map points about MSM point */
00421 float **p_phi, /* coefficients that weight map about MSM points */
00422 int *p_max_phi /* size of phi memory allocation */
00423 ) {
00424 float *phi = NULL;
00425 const int nu = INTERP_PARAMS[msm->interp].nu;
00426 float delta_h, h_delta, t;
00427 int cycle, rmap, diam, nphi;
00428 int i, k;
00429 
00430 *p_h_delta = h_delta = h / delta;
00431 *p_cycle = cycle = n / gcd(m, n);
00432 *p_rmap = rmap = (int) ceilf(h_delta * (nu+1));
00433 
00434 delta_h = delta / h;
00435 diam = 2*rmap + 1;
00436 nphi = diam * cycle;
00437 
00438 if (*p_max_phi < nphi) { /* allocate more memory if we need it */
00439 phi = (float *) realloc(*p_phi, nphi * sizeof(float));
00440 if (NULL == phi) return ERROR(MSMPOT_ERROR_MALLOC);
00441 *p_phi = phi;
00442 *p_max_phi = nphi;
00443 }
00444 ASSERT(*p_phi != NULL);
00445 
00446 for (k = 0; k < cycle; k++) {
00447 float offset = floorf(k * h_delta) * delta_h - k;
00448 phi = *p_phi + k * diam + rmap; /* center of this weight stencil */
00449 switch (msm->interp) {
00450 case MSMPOT_INTERP_CUBIC:
00451 for (i = -rmap; i <= rmap; i++) {
00452 t = fabsf(i * delta_h + offset);
00453 if (t <= 1) {
00454 phi[i] = (1 - t) * (1 + t - 1.5f * t * t);
00455 }
00456 else if (t <= 2) {
00457 phi[i] = 0.5f * (1 - t) * (2 - t) * (2 - t);
00458 }
00459 else {
00460 phi[i] = 0;
00461 }
00462 }
00463 break;
00464 case MSMPOT_INTERP_QUINTIC:
00465 for (i = -rmap; i <= rmap; i++) {
00466 t = fabsf(i * delta_h);
00467 if (t <= 1) {
00468 phi[i] = (1-t*t) * (2-t) * (0.5f + t * (0.25f - (5.f/12)*t));
00469 }
00470 else if (t <= 2) {
00471 phi[i] = (1-t)*(2-t)*(3-t) * ((1.f/6) + t*(0.375f - (5.f/24)*t));
00472 }
00473 else if (t <= 3) {
00474 phi[i] = (1.f/24) * (1-t) * (2-t) * (3-t) * (3-t) * (4-t);
00475 }
00476 else {
00477 phi[i] = 0;
00478 }
00479 }
00480 break;
00481 default:
00482 return ERRMSG(MSMPOT_ERROR_SUPPORT,
00483 "interpolation method not implemented");
00484 } /* end switch on interp */
00485 } /* end loop k over cycles */
00486 return MSMPOT_SUCCESS;
00487 }
00488 
00489 
00490 int setup_domain(Msmpot *msm) {
00491 const float *atom = msm->atom;
00492 const int natoms = msm->natoms;
00493 int n;
00494 float xmin, xmax, ymin, ymax, zmin, zmax;
00495 
00496 /* find extent of atoms */
00497 ASSERT(natoms > 0);
00498 xmin = xmax = atom[ATOM_X(0)];
00499 ymin = ymax = atom[ATOM_Y(0)];
00500 zmin = zmax = atom[ATOM_Z(0)];
00501 for (n = 1; n < natoms; n++) {
00502 float x = atom[ATOM_X(n)];
00503 float y = atom[ATOM_Y(n)];
00504 float z = atom[ATOM_Z(n)];
00505 if (xmin > x) xmin = x;
00506 else if (xmax < x) xmax = x;
00507 if (ymin > y) ymin = y;
00508 else if (ymax < y) ymax = y;
00509 if (zmin > z) zmin = z;
00510 else if (zmax < z) zmax = z;
00511 }
00512 
00513 /* store maximum extent of atoms, regardless of periodicity */
00514 msm->xmin = xmin;
00515 msm->xmax = xmax;
00516 msm->ymin = ymin;
00517 msm->ymax = ymax;
00518 msm->zmin = zmin;
00519 msm->zmax = zmax;
00520 
00521 #if 1
00522 /* domain for non-periodic dimensions is to include both epotmap and atoms */
00523 if ( ! IS_SET_X(msm->isperiodic) ) { /* non-periodic in x */
00524 float lx1 = msm->lx0 + msm->lx; /* contains last epotmap point */
00525 if (xmin >= msm->lx0 && xmax < lx1) {
00526 msm->px = msm->lx; /* assignment can enable factored interpolation */
00527 msm->px0 = msm->lx0;
00528 }
00529 else {
00530 if (xmin > msm->lx0) xmin = msm->lx0;
00531 msm->px0 = xmin;
00532 if (xmax < lx1) {
00533 xmax = lx1;
00534 msm->px = xmax - xmin;
00535 }
00536 else {
00537 msm->px = xmax - xmin + msm->dx;
00538 }
00539 }
00540 }
00541 if ( ! IS_SET_Y(msm->isperiodic) ) { /* non-periodic in y */
00542 float ly1 = msm->ly0 + msm->ly; /* contains last epotmap point */
00543 if (ymin >= msm->ly0 && ymax < ly1) {
00544 msm->py = msm->ly; /* assignment can enable factored interpolation */
00545 msm->py0 = msm->ly0;
00546 }
00547 else {
00548 if (ymin > msm->ly0) ymin = msm->ly0;
00549 msm->py0 = ymin;
00550 if (ymax < ly1) {
00551 ymax = ly1;
00552 msm->py = ymax - ymin;
00553 }
00554 else {
00555 msm->py = ymax - ymin + msm->dy;
00556 }
00557 }
00558 }
00559 if ( ! IS_SET_Z(msm->isperiodic) ) { /* non-periodic in z */
00560 float lz1 = msm->lz0 + msm->lz; /* contains last epotmap point */
00561 if (zmin >= msm->lz0 && zmax < lz1) {
00562 msm->pz = msm->lz; /* assignment can enable factored interpolation */
00563 msm->pz0 = msm->lz0;
00564 }
00565 else {
00566 if (zmin > msm->lz0) zmin = msm->lz0;
00567 msm->pz0 = zmin;
00568 if (zmax < lz1) {
00569 zmax = lz1;
00570 msm->pz = zmax - zmin;
00571 }
00572 else {
00573 msm->pz = zmax - zmin + msm->dz;
00574 }
00575 }
00576 }
00577 #else
00578 /* domain for non-periodic dimensions is to include both epotmap and atoms */
00579 if ( ! IS_SET_X(msm->isperiodic) ) { /* non-periodic in x */
00580 float lx1 = msm->lx0 + msm->lx; /* contains last epotmap point */
00581 if (xmin >= msm->lx0 && xmax < lx1) {
00582 msm->px = msm->lx; /* assignment can enable factored interpolation */
00583 msm->px0 = msm->lx0;
00584 }
00585 else {
00586 if (xmin > msm->lx0) xmin = msm->lx0;
00587 if (xmax < lx1) xmax = lx1;
00588 msm->px = xmax - xmin;
00589 msm->px0 = xmin;
00590 }
00591 }
00592 if ( ! IS_SET_Y(msm->isperiodic) ) { /* non-periodic in y */
00593 float ly1 = msm->ly0 + msm->ly; /* contains last epotmap point */
00594 if (ymin >= msm->ly0 && ymax < ly1) {
00595 msm->py = msm->ly; /* assignment can enable factored interpolation */
00596 msm->py0 = msm->ly0;
00597 }
00598 else {
00599 if (ymin > msm->ly0) ymin = msm->ly0;
00600 if (ymax < ly1) ymax = ly1;
00601 msm->py = ymax - ymin;
00602 msm->py0 = ymin;
00603 }
00604 }
00605 if ( ! IS_SET_Z(msm->isperiodic) ) { /* non-periodic in z */
00606 float lz1 = msm->lz0 + msm->lz; /* contains last epotmap point */
00607 if (zmin >= msm->lz0 && zmax < lz1) {
00608 msm->pz = msm->lz; /* assignment can enable factored interpolation */
00609 msm->pz0 = msm->lz0;
00610 }
00611 else {
00612 if (zmin > msm->lz0) zmin = msm->lz0;
00613 if (zmax < lz1) zmax = lz1;
00614 msm->pz = zmax - zmin;
00615 msm->pz0 = zmin;
00616 }
00617 }
00618 #endif
00619 
00620 return MSMPOT_SUCCESS;
00621 }
00622 
00623 
00624 /*
00625 * Bins are set up based on domain side lengths.
00626 * This is done independently from selecting origin for periodic systems.
00627 */
00628 int setup_bins(Msmpot *msm) {
00629 /* vol is measure of local volume */
00630 float vol = msm->binfill * msm->bindepth / msm->density;
00631 float blen = powf(vol, 1.f/3); /* ideal bin side length */
00632 int maxbin;
00633 int nbx = (int) ceilf(msm->px / blen); /* using ceilf to count bins */
00634 int nby = (int) ceilf(msm->py / blen); /* makes bin vol <= desired vol */
00635 int nbz = (int) ceilf(msm->pz / blen);
00636 
00637 ASSERT(nbx > 0);
00638 ASSERT(nby > 0);
00639 ASSERT(nbz > 0);
00640 
00641 maxbin = nbx * nby * nbz;
00642 if (msm->maxbin < maxbin) { /* grab more memory if we need it */
00643 void *v;
00644 size_t floatsperbin = msm->bindepth * ATOM_SIZE;
00645 v = realloc(msm->bin, maxbin * floatsperbin * sizeof(float));
00646 if (NULL == v) return ERROR(MSMPOT_ERROR_MALLOC);
00647 msm->bin = (float *) v;
00648 v = realloc(msm->bincount, maxbin * sizeof(int));
00649 if (NULL == v) return ERROR(MSMPOT_ERROR_MALLOC);
00650 msm->bincount = (int *) v;
00651 msm->maxbin = maxbin;
00652 
00653 /* for cursor linked list implementation */
00654 v = realloc(msm->first_atom_index, maxbin * sizeof(int));
00655 if (NULL == v) return ERROR(MSMPOT_ERROR_MALLOC);
00656 msm->first_atom_index = (int *) v;
00657 /* */
00658 }
00659 
00660 /* for cursor linked list implmentation */
00661 if (msm->maxatoms < msm->natoms) {
00662 void *v;
00663 v = realloc(msm->next_atom_index, msm->natoms * sizeof(int));
00664 if (NULL == v) return ERROR(MSMPOT_ERROR_MALLOC);
00665 msm->next_atom_index = (int *) v;
00666 msm->maxatoms = msm->natoms;
00667 }
00668 /* */
00669 
00670 msm->nbx = nbx;
00671 msm->nby = nby;
00672 msm->nbz = nbz;
00673 
00674 msm->bx = msm->px / nbx;
00675 msm->by = msm->py / nby;
00676 msm->bz = msm->pz / nbz;
00677 
00678 msm->invbx = 1 / msm->bx;
00679 msm->invby = 1 / msm->by;
00680 msm->invbz = 1 / msm->bz;
00681 
00682 if (msm->maxover < DEFAULT_OVER) {
00683 void *vover = realloc(msm->over, DEFAULT_OVER * ATOM_SIZE * sizeof(float));
00684 if (NULL == vover) return ERROR(MSMPOT_ERROR_MALLOC);
00685 msm->over = (float *) vover;
00686 msm->maxover = DEFAULT_OVER;
00687 }
00688 
00689 return MSMPOT_SUCCESS;
00690 }
00691 
00692 
00693 int setup_origin(Msmpot *msm) {
00694 /* we get to choose the origin (lower-left corner) of the atom domain,
00695 * which can be advantageous to reduce wrapping for periodic boundaries */
00696 
00697 msm->isbinwrap = 0; /* reset flags */
00698 msm->islongcutoff = 0;
00699 
00700 if (IS_SET_X(msm->isperiodic)) {
00701 /* how many bins does cutoff extend? */
00702 int mbx = (int) ceilf(msm->a * msm->invbx);
00703 if (msm->lx < (msm->px - 2*mbx*msm->bx)) {
00704 /* epotmap fits inside of domain with thick shell of bins */
00705 msm->px0 = msm->lx0 - mbx * msm->bx;
00706 }
00707 else {
00708 /* we will have to wrap bin neighborhood */
00709 msm->px0 = msm->lx0;
00710 SET_X(msm->isbinwrap);
00711 if (mbx > msm->nbx) SET_X(msm->islongcutoff); /* wraps more than once */
00712 }
00713 }
00714 
00715 if (IS_SET_Y(msm->isperiodic)) {
00716 /* how many bins does cutoff extend? */
00717 int mby = (int) ceilf(msm->a * msm->invby);
00718 if (msm->ly < (msm->py - 2*mby*msm->by)) {
00719 /* epotmap fits inside of domain with thick shell of bins */
00720 msm->py0 = msm->ly0 - mby * msm->by;
00721 }
00722 else {
00723 /* we will have to wrap bin neighborhood */
00724 msm->py0 = msm->ly0;
00725 SET_Y(msm->isbinwrap);
00726 if (mby > msm->nby) SET_Y(msm->islongcutoff); /* wraps more than once */
00727 }
00728 }
00729 
00730 if (IS_SET_Z(msm->isperiodic)) {
00731 /* how many bins does cutoff extend? */
00732 int mbz = (int) ceilf(msm->a * msm->invbz);
00733 if (msm->lz < (msm->pz - 2*mbz*msm->bz)) {
00734 /* epotmap fits inside of domain with thick shell of bins */
00735 msm->pz0 = msm->lz0 - mbz * msm->bz;
00736 }
00737 else {
00738 /* we will have to wrap bin neighborhood */
00739 msm->pz0 = msm->lz0;
00740 SET_Z(msm->isbinwrap);
00741 if (mbz > msm->nbz) SET_Z(msm->islongcutoff); /* wraps more than once */
00742 }
00743 }
00744 
00745 return MSMPOT_SUCCESS;
00746 }
00747 
00748 
00749 int setup_nonperiodic_hlevelparams_1d(Msmpot *msm,
00750 float len, /* measure to furthest point interpolated from grid */
00751 float *hh, /* determine h */
00752 int *nn, /* determine number grid spacings covering domain */
00753 int *aindex, /* determine smallest lattice index */
00754 int *bindex /* determine largest lattice index */
00755 ) {
00756 const float hmin = msm->hmin; /* minimum bound on h */
00757 const int nu = INTERP_PARAMS[msm->interp].nu; /* interp stencil radius */
00758 /* make sure RH grid point is beyond farthest atom or epotmap point */
00759 int n = (int) floorf(len / hmin) + 1;
00760 *hh = hmin;
00761 *nn = n;
00762 *aindex = -nu;
00763 *bindex = n + nu;
00764 return MSMPOT_SUCCESS;
00765 }
00766 
00767 
00768 int setup_periodic_hlevelparams_1d(Msmpot *msm,
00769 float len, /* domain length */
00770 float *hh, /* determine h */
00771 int *nn, /* determine number grid spacings covering domain */
00772 int *aindex, /* determine smallest lattice index */
00773 int *bindex /* determine largest lattice index */
00774 ) {
00775 const float hmin = msm->hmin; /* minimum bound on h */
00776 const float hmax = 1.5f * hmin;
00777 float h = len;
00778 int n = 1; /* start with one grid point across domain */
00779 while (h >= hmax) {
00780 h *= 0.5f; /* halve h */
00781 n <<= 1; /* double grid points */
00782 }
00783 if (h < hmin) {
00784 if (n < 4) { /* error: either len is too small or hmin is too large */
00785 return ERRMSG(MSMPOT_ERROR_PARAM,
00786 "ratio of domain length to hmin is too small");
00787 }
00788 h *= (4.f/3); /* scale h by 4/3 */
00789 n >>= 2; /* scale n by 3/4 */
00790 n *= 3;
00791 }
00792 /* now we have: hmin <= h < hmax */
00793 /* now we have: n is power of two times no more than one power of 3 */
00794 *hh = h;
00795 *nn = n;
00796 *aindex = 0;
00797 *bindex = n-1;
00798 return MSMPOT_SUCCESS;
00799 }
00800 
00801 
00802 int setup_hierarchy(Msmpot *msm) {
00803 const int nu = INTERP_PARAMS[msm->interp].nu;
00804 const int omega = INTERP_PARAMS[msm->interp].omega;
00805 const int split = msm->split;
00806 const int ispx = IS_SET_X(msm->isperiodic);
00807 const int ispy = IS_SET_Y(msm->isperiodic);
00808 const int ispz = IS_SET_Z(msm->isperiodic);
00809 const int ispany = IS_SET_ANY(msm->isperiodic);
00810 const float a = msm->a;
00811 float hx, hy, hz;
00812 float scaling;
00813 
00814 floatGrid *p = NULL;
00815 int ia, ib, ja, jb, ka, kb, ni, nj, nk;
00816 int nx, ny, nz; /* counts the grid points that span just the domain */
00817 
00818 int i, j, k, n;
00819 int index;
00820 int level, toplevel, nlevels, maxlevels;
00821 int lastnelems = 1;
00822 int isclamped = 0;
00823 int done, alldone;
00824 
00825 int err = 0;
00826 
00827 if (ispx) {
00828 err = setup_periodic_hlevelparams_1d(msm, msm->px, &hx, &nx, &ia, &ib);
00829 }
00830 else {
00831 float xmax = msm->lx0 + msm->lx - msm->dx; /* furthest epotmap point */
00832 if (xmax < msm->xmax) xmax = msm->xmax; /* furthest atom */
00833 err = setup_nonperiodic_hlevelparams_1d(msm, xmax - msm->px0,
00834 &hx, &nx, &ia, &ib);
00835 }
00836 if (err) return ERROR(err);
00837 
00838 if (ispy) {
00839 err = setup_periodic_hlevelparams_1d(msm, msm->py, &hy, &ny, &ja, &jb);
00840 }
00841 else {
00842 float ymax = msm->ly0 + msm->ly - msm->dy; /* furthest epotmap point */
00843 if (ymax < msm->ymax) ymax = msm->ymax; /* furthest atom */
00844 err = setup_nonperiodic_hlevelparams_1d(msm, ymax - msm->py0,
00845 &hy, &ny, &ja, &jb);
00846 }
00847 if (err) return ERROR(err);
00848 
00849 if (ispz) {
00850 err = setup_periodic_hlevelparams_1d(msm, msm->pz, &hz, &nz, &ka, &kb);
00851 }
00852 else {
00853 float zmax = msm->lz0 + msm->lz - msm->dz; /* furthest epotmap point */
00854 if (zmax < msm->zmax) zmax = msm->zmax; /* furthest atom */
00855 err = setup_nonperiodic_hlevelparams_1d(msm, zmax - msm->pz0,
00856 &hz, &nz, &ka, &kb);
00857 }
00858 if (err) return ERROR(err);
00859 
00860 msm->hx = hx;
00861 msm->hy = hy;
00862 msm->hz = hz;
00863 
00864 msm->nx = nx;
00865 msm->ny = ny;
00866 msm->nz = nz;
00867 
00868 ni = ib - ia + 1;
00869 nj = jb - ja + 1;
00870 nk = kb - ka + 1;
00871 
00872 /* allocate temp buffer space for factored grid transfer */
00873 n = (nk > omega ? nk : omega); /* row along z-dimension */
00874 if (msm->max_lzd < n) {
00875 float *t;
00876 t = (float *) realloc(msm->lzd, n * sizeof(float));
00877 if (NULL == t) return ERROR(MSMPOT_ERROR_MALLOC);
00878 msm->lzd = t;
00879 msm->max_lzd = n;
00880 }
00881 n *= (nj > omega ? nj : omega); /* plane along yz-dimensions */
00882 if (msm->max_lyzd < n) {
00883 float *t;
00884 t = (float *) realloc(msm->lyzd, n * sizeof(float));
00885 if (NULL == t) return ERROR(MSMPOT_ERROR_MALLOC);
00886 msm->lyzd = t;
00887 msm->max_lyzd = n;
00888 }
00889 
00890 nlevels = msm->nlevels;
00891 if (nlevels <= 0) {
00892 /* automatically set number of levels */
00893 n = ni;
00894 if (n < nj) n = nj;
00895 if (n < nk) n = nk;
00896 for (maxlevels = 1; n > 0; n >>= 1) maxlevels++;
00897 nlevels = maxlevels;
00898 if ( ! ispany ) { /* no periodicity */
00899 int omega3 = omega * omega * omega;
00900 int nhalf = (int) sqrtf((float) ni*nj*nk); /* scale down for performance? */
00901 lastnelems = (nhalf > omega3 ? nhalf : omega3);
00902 isclamped = 1;
00903 }
00904 }
00905 else {
00906 /* user-defined number of levels */
00907 maxlevels = nlevels;
00908 }
00909 
00910 /* allocate any additional levels that may be needed */
00911 if (msm->maxlevels < maxlevels) {
00912 void *vqh, *veh, *vgc;
00913 vqh = realloc(msm->qh, maxlevels * sizeof(floatGrid));
00914 if (NULL == vqh) return ERROR(MSMPOT_ERROR_MALLOC);
00915 veh = realloc(msm->eh, maxlevels * sizeof(floatGrid));
00916 if (NULL == veh) return ERROR(MSMPOT_ERROR_MALLOC);
00917 vgc = realloc(msm->gc, maxlevels * sizeof(floatGrid));
00918 if (NULL == vgc) return ERROR(MSMPOT_ERROR_MALLOC);
00919 msm->qh = (floatGrid *) vqh;
00920 msm->eh = (floatGrid *) veh;
00921 msm->gc = (floatGrid *) vgc;
00922 /* initialize the newest grids appended to array */
00923 for (level = msm->maxlevels; level < maxlevels; level++) {
00924 GRID_INIT( &(msm->qh[level]) );
00925 GRID_INIT( &(msm->eh[level]) );
00926 GRID_INIT( &(msm->gc[level]) );
00927 }
00928 msm->maxlevels = maxlevels;
00929 }
00930 
00931 level = 0;
00932 done = 0;
00933 alldone = 0;
00934 do {
00935 GRID_RESIZE( &(msm->qh[level]), ia, ni, ja, nj, ka, nk);
00936 GRID_RESIZE( &(msm->eh[level]), ia, ni, ja, nj, ka, nk);
00937 
00938 if (++level == nlevels) done |= 0x07; /* user limit on levels */
00939 
00940 alldone = (done == 0x07); /* make sure all dimensions are done */
00941 
00942 if (isclamped) {
00943 int nelems = ni * nj * nk;
00944 if (nelems <= lastnelems) done |= 0x07;
00945 }
00946 
00947 if (ispx) {
00948 ni >>= 1;
00949 ib = ni-1;
00950 if (ni & 1) done |= 0x07; /* == 3 or 1 */
00951 else if (ni == 2) done |= 0x01; /* can do one more */
00952 }
00953 else {
00954 ia = -((-ia+1)/2) - nu;
00955 ib = (ib+1)/2 + nu;
00956 ni = ib - ia + 1;
00957 if (ni <= omega) done |= 0x01; /* can do more restrictions */
00958 }
00959 
00960 if (ispy) {
00961 nj >>= 1;
00962 jb = nj-1;
00963 if (nj & 1) done |= 0x07; /* == 3 or 1 */
00964 else if (nj == 2) done |= 0x02; /* can do one more */
00965 }
00966 else {
00967 ja = -((-ja+1)/2) - nu;
00968 jb = (jb+1)/2 + nu;
00969 nj = jb - ja + 1;
00970 if (nj <= omega) done |= 0x02; /* can do more restrictions */
00971 }
00972 
00973 if (ispz) {
00974 nk >>= 1;
00975 kb = nk-1;
00976 if (nk & 1) done |= 0x07; /* == 3 or 1 */
00977 else if (nk == 2) done |= 0x04; /* can do one more */
00978 }
00979 else {
00980 ka = -((-ka+1)/2) - nu;
00981 kb = (kb+1)/2 + nu;
00982 nk = kb - ka + 1;
00983 if (nk <= omega) done |= 0x04; /* can do more restrictions */
00984 }
00985 
00986 } while ( ! alldone );
00987 msm->nlevels = level;
00988 
00989 toplevel = (ispany ? msm->nlevels : msm->nlevels - 1);
00990 
00991 /* ellipsoid axes for lattice cutoff weights */
00992 ni = (int) ceilf(2*a/hx) - 1;
00993 nj = (int) ceilf(2*a/hy) - 1;
00994 nk = (int) ceilf(2*a/hz) - 1;
00995 scaling = 1;
00996 for (level = 0; level < toplevel; level++) {
00997 p = &(msm->gc[level]);
00998 GRID_RESIZE(p, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
00999 
01000 if (0 == level) {
01001 index = 0;
01002 for (k = -nk; k <= nk; k++) {
01003 for (j = -nj; j <= nj; j++) {
01004 for (i = -ni; i <= ni; i++) {
01005 float s, t, gs, gt, g;
01006 s = ( (i*hx)*(i*hx) + (j*hy)*(j*hy) + (k*hz)*(k*hz) ) / (a*a);
01007 t = 0.25f * s;
01008 if (t >= 1) {
01009 g = 0;
01010 }
01011 else if (s >= 1) {
01012 gs = 1/sqrtf(s);
01013 SPOLY(&gt, t, split);
01014 g = (gs - 0.5f * gt) / a;
01015 }
01016 else {
01017 SPOLY(&gs, s, split);
01018 SPOLY(&gt, t, split);
01019 g = (gs - 0.5f * gt) / a;
01020 }
01021 GRID_INDEX_CHECK(p, i, j, k);
01022 ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
01023 p->buffer[index] = g;
01024 index++;
01025 }
01026 }
01027 } /* end loops over k-j-i */
01028 }
01029 else {
01030 /* set each level as scaling of h-level */
01031 const floatGrid *first = &(msm->gc[0]);
01032 scaling *= 0.5f;
01033 index = 0;
01034 for (k = -nk; k <= nk; k++) {
01035 for (j = -nj; j <= nj; j++) {
01036 for (i = -ni; i <= ni; i++) {
01037 GRID_INDEX_CHECK(p, i, j, k);
01038 ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
01039 p->buffer[index] = scaling * first->buffer[index];
01040 index++;
01041 }
01042 }
01043 }
01044 }
01045 } /* end loop over levels */
01046 
01047 if (toplevel < msm->nlevels) {
01048 /* nonperiodic in all dimensions,
01049 * calculate top level weights, ellipsoid axes are length of lattice */
01050 const floatGrid *qhtop = &(msm->qh[toplevel]);
01051 ni = qhtop->ni - 1;
01052 nj = qhtop->nj - 1;
01053 nk = qhtop->nk - 1;
01054 p = &(msm->gc[toplevel]);
01055 GRID_RESIZE(p, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
01056 scaling *= 0.5f;
01057 index = 0;
01058 for (k = -nk; k <= nk; k++) {
01059 for (j = -nj; j <= nj; j++) {
01060 for (i = -ni; i <= ni; i++) {
01061 float s, gs;
01062 s = ( (i*hx)*(i*hx) + (j*hy)*(j*hy) + (k*hz)*(k*hz) ) / (a*a);
01063 if (s >= 1) {
01064 gs = 1/sqrtf(s);
01065 }
01066 else {
01067 SPOLY(&gs, s, split);
01068 }
01069 GRID_INDEX_CHECK(p, i, j, k);
01070 ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
01071 p->buffer[index] = scaling * gs/a;
01072 index++;
01073 }
01074 }
01075 } /* end loops over k-j-i for coarsest level weights */
01076 }
01077 
01078 return 0;
01079 }