/* * fnonbond.c * * Routines to evaluate cutoff nonbonded forces. */ #include #include #include #include #include #include "mdutil.h" #include "fnonbond.h" /* * typedefs */ typedef void (*FunptrComputePair)(Force *, MD_Int, MD_Int, double, MD_Dvec *, MD_Dvec *, MD_Dvec *); /* * prototypes for internal functions */ /* for initialization */ static MD_Errcode init_atombox(Force *force); /* different computational strategies for nonbonded forces */ /* depends on electrostatic evaluation policy */ static MD_Errcode compute_cutoff(Force *force); #if 0 static MD_Errcode compute_mg(Force *force); static MD_Errcode compute_mgrid(Force *force); static MD_Errcode compute_dpmta(Force *force); #endif static MD_Errcode compute_full_direct(Force *force); /* geometric hashing of atoms */ /* returns error code if system has expanded to unreasonable size */ static MD_Errcode hash_atoms(Force *force); /* subtract out elec interactions that should have been excluded */ static void subtract_excl(Force *force); /* compute pairwise interactions between adjacent atom boxes */ /* give appropriate compute_XXpair function */ static void compute_atombox(Force *force, FunptrComputePair compute_pair); /* compute pairwise interaction between atoms */ /* nbpair = both vdw and elec forces */ static void compute_nbpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j); static void compute_vdwpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j); static void compute_elecpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j); static void subtract_elecpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j); static void subtract_elec14pair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j); /* * macros */ /* need to have nxbox and nybox be local (and set correctly) */ #undef BOX_INDEX #define BOX_INDEX(x,y,z) ((x) + nxbox*((y) + nybox*(z))) /* * external functions */ MD_Errcode deven_nonbonded_init(Force *force) { Data *data = force->data; void *tmp; double dis, c; MD_Int k; #if 0 if (data->elec_eval_policy == ELEC_EVAL_MG && data->is_sphericalbc != ON) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_simparam, "deven_nonbonded_init: spherical boundary conditions required for MG"); return MD_FAIL; } if (data->elec_eval_policy == ELEC_EVAL_DPMTA && data->is_sphericalbc != ON) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_simparam, "deven_nonbonded_init: spherical boundary conditions required for DPMTA"); return MD_FAIL; } #endif switch (data->elec_eval_policy) { case ELEC_EVAL_CUTOFF: force->compute_nonbonded = compute_cutoff; break; case ELEC_EVAL_FULLDIRECT: force->compute_nonbonded = compute_full_direct; break; #if 0 case ELEC_EVAL_MG: /* validate mg simparams */ dis = ((data->sphericalbc_r1 > data->sphericalbc_r2) ? data->sphericalbc_r1 : data->sphericalbc_r2) * 2.0; if (data->mg_divisions * data->mg_gridsize < dis) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_simparam, "deven_nonbonded_init: mg_divisions * mg_gridsize must be " "greater than diameter of bounding sphere"); return MD_FAIL; } /* set pointer to nonbonded force computation routine */ force->compute_nonbonded = compute_mg; /* allocate memory for MG position and charge arrays */ /* * This really sucks! Have to allocate an additional position * array because MG routine uses pos[0] to store its "origin". * This will also require an additional copy from our own * position array before each electrostatic evaluation. */ if (force->mg_pos != NULL) freeMatrix(force->mg_pos); if ((force->mg_pos = matrix(data->natoms + 1, 3)) == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to matrix (mg)"); return MD_FAIL; } tmp = realloc(force->mg_q, data->natoms * sizeof(double)); if (tmp == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to realloc"); return MD_FAIL; } force->mg_q = (double *) tmp; /* initialize the array of charges */ /* charge requires a conversion factor */ c = sqrt(force->elec_const / 321.008749); for (k = 0; k < data->natoms; k++) { force->mg_q[k] = c * data->atom[k].q; } /* set MG parameters */ force->mg_params.N = data->natoms; force->mg_params.M = data->mg_divisions; force->mg_params.h = data->mg_gridsize; force->mg_params.Level = data->mg_levels; force->mg_params.a = data->mg_cutoff; /* call MG initialization routine */ /* ("dis" is a useless parameter that gets zeroed) */ if (force->mg_f != NULL) freeMatrix(force->mg_f); mg_setUp(&(force->mg_f), &dis, force->mg_pos, &(force->mg_params)); if (force->mg_f == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to mg_setUp"); return MD_FAIL; } /* override the "center" calculated by mg_setUp */ /* it should instead be the center of our bounding sphere */ force->mg_pos[0][0] = data->sphericalbc_center.x; force->mg_pos[0][1] = data->sphericalbc_center.y; force->mg_pos[0][2] = data->sphericalbc_center.z; break; case ELEC_EVAL_MGRID: /* one of MGgridsize and MGlength is optional */ if (data->mg_divisions == 0) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_simparam, "deven_nonbonded_init: mg_divisions must be greater than 0"); return MD_FAIL; } if (data->mg_gridsize == 0) { data->mg_gridsize = data->mg_length / data->mg_divisions; } else if (data->mg_length == 0) { data->mg_length = data->mg_divisions * data->mg_gridsize; } /* validate mgrid simparams */ dis = ((data->sphericalbc_r1 > data->sphericalbc_r2) ? data->sphericalbc_r1 : data->sphericalbc_r2) * 2.0; if (data->mg_divisions * data->mg_gridsize < dis) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_simparam, "deven_nonbonded_init: mg_divisions * mg_gridsize must be " "greater than diameter of bounding sphere"); return MD_FAIL; } /* set pointer to nonbonded force computation routine */ force->compute_nonbonded = compute_mgrid; /* allocate memory for mgrid charge array */ tmp = realloc(force->mgrid_q, data->natoms * sizeof(double)); if (tmp == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to realloc"); return MD_FAIL; } force->mgrid_q = (double *) tmp; /* initialize the array of charges */ /* charge requires a conversion factor */ c = sqrt(force->elec_const); for (k = 0; k < data->natoms; k++) { force->mgrid_q[k] = c * data->atom[k].q; } /* set mgrid params */ force->mgrid_params.is_periodic = 0; force->mgrid_params.nlevels = data->mg_levels; force->mgrid_params.ndivisions = data->mg_divisions; force->mgrid_params.gridsize = data->mg_gridsize; force->mgrid_params.cutoff = data->mg_cutoff; force->mgrid_params.length = data->mg_length; force->mgrid_params.center.x = data->sphericalbc_center.x; force->mgrid_params.center.y = data->sphericalbc_center.y; force->mgrid_params.center.z = data->sphericalbc_center.z; force->mgrid_params.switchtype = data->mg_switchtype; /* call mgrid creation routine */ force->mgrid = mgrid_create(&(force->mgrid_params)); if (force->mgrid == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_mgrid, "deven_nonbonded_init call to mgrid_create failed"); return MD_FAIL; } break; case ELEC_EVAL_DPMTA: /* set pointer to nonbonded force computation routine */ force->compute_nonbonded = compute_dpmta; /* initialize DPMTA parameters */ dis = ((data->sphericalbc_r1 > data->sphericalbc_r2) ? data->sphericalbc_r1 : data->sphericalbc_r2) * 2.0 + 2.0; force->pmta_init_data.v1.x = dis; force->pmta_init_data.v2.y = dis; force->pmta_init_data.v3.z = dis; force->pmta_init_data.v1.y = 0; force->pmta_init_data.v1.z = 0; force->pmta_init_data.v2.x = 0; force->pmta_init_data.v2.z = 0; force->pmta_init_data.v3.x = 0; force->pmta_init_data.v3.y = 0; force->pmta_init_data.cellctr.x = data->sphericalbc_center.x; force->pmta_init_data.cellctr.y = data->sphericalbc_center.y; force->pmta_init_data.cellctr.z = data->sphericalbc_center.z; force->pmta_init_data.nprocs = 1; force->pmta_init_data.nlevels = data->dpmta_levels; force->pmta_init_data.mp = data->dpmta_mpterms; force->pmta_init_data.mp_lj = 0; force->pmta_init_data.mp_vir = 0; force->pmta_init_data.fft = data->dpmta_is_fft; force->pmta_init_data.fftblock = data->dpmta_blocking_fft; force->pmta_init_data.pbc = 0; force->pmta_init_data.kterm = 0; force->pmta_init_data.theta = data->dpmta_theta; force->pmta_init_data.calling_num = 0; force->pmta_init_data.calling_tids = NULL; force->pmta_init_data.loadstep = 0; force->pmta_init_data.loadweight = 0; if (PMTAinit(&(force->pmta_init_data), NULL)) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_dpmta, "deven_nonbonded_init call to PMTAinit"); return MD_FAIL; } /* allocate memory needed for calls to DPMTA */ tmp = realloc(force->pmta_particle, data->natoms * sizeof(PmtaParticle)); if (tmp == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to realloc"); return MD_FAIL; } force->pmta_particle = (PmtaParticle *) tmp; tmp = realloc(force->pmta_result, data->natoms * sizeof(PmtaPartInfo)); if (tmp == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to realloc"); return MD_FAIL; } force->pmta_result = (PmtaPartInfo *) tmp; /* initialize the charge data */ /* charge requires a conversion factor */ c = sqrt(force->elec_const); for (k = 0; k < data->natoms; k++) { force->pmta_particle[k].q = c * data->atom[k].q; } break; #endif default: MD_ERR(data->sim, MD_engine_name(data->sim), data->err_simparam, "deven_nonbonded_init: illegal electrostatic evaluation policy"); return MD_FAIL; } return init_atombox(force); } void deven_nonbonded_destroy(Force *force) { AtomBox *atombox = force->atombox; const MD_Int numbox = force->numbox; MD_Int i; #if 0 /* clean up MG */ mg_cleanUp(force->mg_f, force->mg_pos, force->mg_q); force->mg_f = force->mg_pos = NULL; force->mg_q = NULL; /* clean up mgrid */ if (force->mgrid) mgrid_destroy(force->mgrid); free(force->mgrid_q); force->mgrid_q = NULL; /* clean up DPMTA */ free(force->pmta_particle); force->pmta_particle = NULL; free(force->pmta_result); force->pmta_result = NULL; #endif /* free memory used by each atom box */ if (atombox != NULL) { for (i = 0; i < numbox; i++) free(atombox[i].atomid); free(atombox); } } /****************************************************************************/ /** **/ /** internal helper functions **/ /** **/ /****************************************************************************/ /* * initialize atom boxes * fails if memory can't be allocated */ MD_Errcode init_atombox(Force *force) { Data *data = force->data; const MD_Int natoms = data->natoms; MD_Dvec *pos = data->pos; double radius, xmin, xmax, ymin, ymax, zmin, zmax; double xrange, yrange, zrange, xboxsize, yboxsize, zboxsize; double inv_xboxsize, inv_yboxsize, inv_zboxsize; double xadjust, yadjust, zadjust; AtomBox *atombox; MD_Int *nbr; MD_Int nxbox, nybox, nzbox, numbox, numnbrs; MD_Int index, i, j, k, i_plus, i_minus, j_plus, j_minus, k_plus; char buf[BUFSIZE]; /* find extent of system */ if (force->is_sphericalbc) { /* use radius of largest sphere for extent */ if (force->r2 > force->r1) radius = force->r2; else radius = force->r1; xmin = force->center.x - radius; xmax = force->center.x + radius; ymin = force->center.y - radius; ymax = force->center.y + radius; zmin = force->center.z - radius; zmax = force->center.z + radius; } else { /* find max and min of all atom positions */ xmin = xmax = pos[0].x; ymin = ymax = pos[0].y; zmin = zmax = pos[0].z; for (i = 1; i < natoms; i++) { if (pos[i].x < xmin) xmin = pos[i].x; else if (pos[i].x > xmax) xmax = pos[i].x; if (pos[i].y < ymin) ymin = pos[i].y; else if (pos[i].y > ymax) ymax = pos[i].y; if (pos[i].z < zmin) zmin = pos[i].z; else if (pos[i].z > zmax) zmax = pos[i].z; } } xrange = xmax - xmin; yrange = ymax - ymin; zrange = zmax - zmin; xboxsize = yboxsize = zboxsize = data->cutoff; inv_xboxsize = 1.0 / xboxsize; inv_yboxsize = 1.0 / yboxsize; inv_zboxsize = 1.0 / zboxsize; /* determine number of boxes in three directions */ /* so boxes entirely encompass extent of system */ nxbox = (MD_Int) (xrange * inv_xboxsize); if (nxbox * xboxsize < xrange || nxbox == 0) nxbox++; nybox = (MD_Int) (yrange * inv_yboxsize); if (nybox * yboxsize < yrange || nybox == 0) nybox++; nzbox = (MD_Int) (zrange * inv_zboxsize); if (nzbox * zboxsize < zrange || nzbox == 0) nzbox++; numbox = nxbox * nybox * nzbox; /* readjust the range based on the boxes */ xadjust = 0.5 * (nxbox * xboxsize - xrange); xmin -= xadjust; xmax += xadjust; yadjust = 0.5 * (nybox * yboxsize - yrange); ymin -= yadjust; ymax += yadjust; zadjust = 0.5 * (nzbox * zboxsize - zrange); zmin -= zadjust; zmax += zadjust; /* allocate the atom boxes */ atombox = (AtomBox *) calloc(numbox, sizeof(AtomBox)); if (atombox == NULL) return MD_ERR_MEMALLOC; for (i = 0; i < numbox; i++) { atombox[i].atomid = (MD_Int *) malloc(MIN_ATOMIDS * sizeof(MD_Int)); if (atombox[i].atomid == NULL) { MD_ERR(data->sim, MD_engine_name(data->sim), MD_ERR_MEMALLOC, "deven_nonbonded_init call to malloc"); return MD_FAIL; } atombox[i].len = MIN_ATOMIDS; } /* init each box */ for (k = 0; k < nzbox; k++) { for (j = 0; j < nybox; j++) { for (i = 0; i < nxbox; i++) { index = BOX_INDEX(i, j, k); /* set box ID and extent */ atombox[index].nx = i; atombox[index].ny = j; atombox[index].nz = k; atombox[index].min.x = xmin + i * xboxsize; atombox[index].max.x = xmin + (i+1) * xboxsize; atombox[index].min.y = ymin + j * yboxsize; atombox[index].max.y = ymin + (j+1) * yboxsize; atombox[index].min.z = zmin + k * zboxsize; atombox[index].max.z = zmin + (k+1) * zboxsize; /* adjust extent of outer boxes */ /* note that 0 == nxbox-1 when only one x box, etc. */ if (i == 0) atombox[index].min.x = -DBL_MAX; if (i == nxbox-1) atombox[index].max.x = DBL_MAX; if (j == 0) atombox[index].min.y = -DBL_MAX; if (j == nybox-1) atombox[index].max.y = DBL_MAX; if (k == 0) atombox[index].min.z = -DBL_MAX; if (k == nzbox-1) atombox[index].max.z = DBL_MAX; /* create neighbor lists for nonperiodic system */ nbr = atombox[index].nbr; numnbrs = 0; i_plus = (i + 1 < nxbox); i_minus = (i > 0); j_plus = (j + 1 < nybox); j_minus = (j > 0); k_plus = (k + 1 < nzbox); if (i_plus) nbr[numnbrs++] = BOX_INDEX(i+1, j, k); if (j_plus) { if (i_minus) nbr[numnbrs++] = BOX_INDEX(i-1, j+1, k); nbr[numnbrs++] = BOX_INDEX(i, j+1, k); if (i_plus) nbr[numnbrs++] = BOX_INDEX(i+1, j+1, k); } if (k_plus) { if (j_minus) { if (i_minus) nbr[numnbrs++] = BOX_INDEX(i-1, j-1, k+1); nbr[numnbrs++] = BOX_INDEX(i, j-1, k+1); if (i_plus) nbr[numnbrs++] = BOX_INDEX(i+1, j-1, k+1); } if (i_minus) nbr[numnbrs++] = BOX_INDEX(i-1, j, k+1); nbr[numnbrs++] = BOX_INDEX(i, j, k+1); if (i_plus) nbr[numnbrs++] = BOX_INDEX(i+1, j, k+1); if (j_plus) { if (i_minus) nbr[numnbrs++] = BOX_INDEX(i-1, j+1, k+1); nbr[numnbrs++] = BOX_INDEX(i, j+1, k+1); if (i_plus) nbr[numnbrs++] = BOX_INDEX(i+1, j+1, k+1); } } atombox[index].numnbrs = numnbrs; } } } force->atombox = atombox; force->numbox = numbox; force->nxbox = nxbox; force->nybox = nybox; force->nzbox = nzbox; force->xmin = xmin; force->ymin = ymin; force->zmin = zmin; force->inv_xboxsize = inv_xboxsize; force->inv_yboxsize = inv_yboxsize; force->inv_zboxsize = inv_zboxsize; force->xextent = xrange; force->yextent = yrange; force->zextent = zrange; return 0; } /* * simple nonbonded cutoff algorithm: * always hash before doing a box force evaluation */ MD_Errcode compute_cutoff(Force *force) { if (hash_atoms(force)) return MD_FAIL; compute_atombox(force, compute_nbpair); return 0; } /* * use this for debugging purposes * still need to hash atoms for cutoff VDW evaluation */ MD_Errcode compute_full_direct(Force *force) { MD_Dvec r_ij; MD_Dvec *pos = force->data->pos; MD_Dvec *f = force->data->f; double rlen2; MD_Int natoms = force->data->natoms; MD_Int i, j; /* compute elec forces and potential between all pairs of atoms */ for (j = 1; j < natoms; j++) { for (i = 0; i < j; i++) { r_ij.x = pos[j].x - pos[i].x; r_ij.y = pos[j].y - pos[i].y; r_ij.z = pos[j].z - pos[i].z; rlen2 = r_ij.x * r_ij.x + r_ij.y * r_ij.y + r_ij.z * r_ij.z; compute_elecpair(force, i, j, rlen2, &r_ij, &(f[i]), &(f[j])); } } /* subtract force (and potential) elec contributions from excluded pairs */ subtract_excl(force); /* hash atoms to compute cutoff VDW force */ if (hash_atoms(force)) return MD_FAIL; compute_atombox(force, compute_vdwpair); return 0; } #if 0 /* * call Ismail Tezcan's MG N-body solver library * subtract out electrostatic exclusions * use atom hashing for cutoff VDW evaluation */ MD_Errcode compute_mg(Force *force) { Data *data = force->data; MD_Int natoms = data->natoms; MD_Int k; /* zero MG force vector array */ memset(force->mg_f[0], 0, natoms * sizeof(MD_Dvec)); /* copy our current positions */ memcpy(force->mg_pos[1], data->pos, natoms * sizeof(MD_Dvec)); /* call MG n-body solver */ mg_compute(force->mg_f, &(force->pe_elec), force->mg_pos, force->mg_q, &(force->mg_params)); /* accumulate elec forces into total force, convert to internal units */ for (k = 0; k < natoms; k++) { data->f[k].x += force->mg_f[k][0]; data->f[k].y += force->mg_f[k][1]; data->f[k].z += force->mg_f[k][2]; } /* subtract force (and potential) elec contributions from excluded pairs */ subtract_excl(force); /* hash atoms to compute cutoff VDW force */ if (hash_atoms(force)) return MD_FAIL; compute_atombox(force, compute_vdwpair); return 0; } /* * call my mgrid N-body solver library * subtract out electrostatic exclusions * use atom hashing for cutoff VDW evaluation */ MD_Errcode compute_mgrid(Force *force) { Data *data = force->data; /* set system params to mgrid */ force->mgrid_system.force = (MgridVector *) data->f; force->mgrid_system.position = (MgridVector *) data->pos; force->mgrid_system.charge = force->mgrid_q; force->mgrid_system.natoms = data->natoms; /* call mgrid N-body solver */ if (mgrid_compute(force->mgrid, &(force->mgrid_system))) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_mgrid, "compute_mgrid call to mgrid_compute"); return MD_FAIL; } force->pe_elec = force->mgrid_system.potential; /* subtract force (and potential) elec contributions from excluded pairs */ subtract_excl(force); /* hash atoms to compute cutoff VDW force */ if (hash_atoms(force)) return MD_FAIL; compute_atombox(force, compute_vdwpair); return 0; } /* * call DPMTA * subtract out electrostatic exclusions * use atom hashing for cutoff VDW evaluation */ MD_Errcode compute_dpmta(Force *force) { Data *data = force->data; MD_Int natoms = data->natoms; MD_Int k; double pe_elec = 0; /* zero DPMTA force and elec potential array */ memset(force->pmta_result, 0, natoms * sizeof(PmtaPartInfo)); /* copy our current positions */ for (k = 0; k < natoms; k++) { force->pmta_particle[k].p.x = data->pos[k].x; force->pmta_particle[k].p.y = data->pos[k].y; force->pmta_particle[k].p.z = data->pos[k].z; } /* call DPMTA */ if (PMTAforce(natoms, force->pmta_particle, force->pmta_result, NULL)) { MD_ERR(data->sim, MD_engine_name(data->sim), data->err_dpmta, "compute_dpmta call to PMTAforce"); return MD_FAIL; } /* accumulate forces and total elec potential from individual particles */ for (k = 0; k < natoms; k++) { data->f[k].x += force->pmta_result[k].f.x; data->f[k].y += force->pmta_result[k].f.y; data->f[k].z += force->pmta_result[k].f.z; pe_elec += force->pmta_result[k].v; } /* counted potentials twice, have to divide by 2 */ force->pe_elec = 0.5 * pe_elec; /* subtract force (and potential) elec contributions from excluded pairs */ subtract_excl(force); /* hash atoms to compute cutoff VDW force */ if (hash_atoms(force)) return MD_FAIL; compute_atombox(force, compute_vdwpair); return 0; } #endif /* * perform geometric hashing on atoms, place into boxes * (assume nonperiodic boundaries) * fails if realloc of arrays fail or if system size expands unreasonably */ MD_Errcode hash_atoms(Force *force) { MD_Dvec *pos = force->data->pos; const MD_Int natoms = force->data->natoms; const MD_Int numbox = force->numbox; const MD_Int nxbox = force->nxbox; const MD_Int nybox = force->nybox; const MD_Int nzbox = force->nzbox; double xmin = force->xmin; double ymin = force->ymin; double zmin = force->zmin; double xmax, ymax, zmax; const double inv_xboxsize = force->inv_xboxsize; const double inv_yboxsize = force->inv_yboxsize; const double inv_zboxsize = force->inv_zboxsize; AtomBox *box = force->atombox; MD_Int i, k, kx, ky, kz; /* recompute system min if not using spherical boundaries */ /* do this in case system tries to drift */ /* if system becomes 10 times its normal size, return an error message */ if (!force->is_sphericalbc) { xmin = xmax = pos[0].x; ymin = ymax = pos[0].y; zmin = zmax = pos[0].z; for (i = 1; i < natoms; i++) { if (pos[i].x < xmin) xmin = pos[i].x; else if (pos[i].x > xmax) xmax = pos[i].x; if (pos[i].y < ymin) ymin = pos[i].y; else if (pos[i].y > ymax) ymax = pos[i].y; if (pos[i].z < zmin) zmin = pos[i].z; else if (pos[i].z > zmax) zmax = pos[i].z; } if (xmax - xmin > 10 * force->xextent || ymax - ymin > 10 * force->yextent || zmax - zmin > 10 * force->zextent) { MD_ERR(force->data->sim, MD_engine_name(force->data->sim), force->data->err_unstable, "hash_atoms found system size expanded 10 times original extent"); return MD_FAIL; } force->xmin = xmin; force->ymin = ymin; force->zmin = zmin; } /* init each atom box */ for (k = 0; k < numbox; k++) { /* check to see if number of atoms warrants reducing array size */ if (box[k].len > MIN_ATOMIDS && box[k].num <= box[k].len / 4) { void *tmp = realloc(box[k].atomid, sizeof(MD_Int) * box[k].len / 2); if (tmp == NULL) { MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "hash_atoms call to realloc"); return MD_FAIL; } box[k].atomid = (MD_Int *) tmp; box[k].len /= 2; } /* clear number of atoms in each box for rehashing */ box[k].num = 0; } /* place each atom into its correct box */ for (i = 0; i < natoms; i++) { /* hash the atom */ kx = (int)((pos[i].x - xmin) * inv_xboxsize); if (kx < 0) kx = 0; else if (kx >= nxbox) kx = nxbox-1; ky = (int)((pos[i].y - ymin) * inv_yboxsize); if (ky < 0) ky = 0; else if (ky >= nybox) ky = nybox-1; kz = (int)((pos[i].z - zmin) * inv_zboxsize); if (kz < 0) kz = 0; else if (kz >= nzbox) kz = nzbox-1; k = BOX_INDEX(kx, ky, kz); /* resize atomid array for this box if needed */ if (box[k].num == box[k].len) { void *tmp = realloc(box[k].atomid, sizeof(MD_Int) * box[k].len * 2); if (tmp == NULL) { MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "hash_atoms call to realloc"); return MD_FAIL; } box[k].atomid = (MD_Int *) tmp; box[k].len *= 2; } box[k].atomid[ box[k].num++ ] = i; } return 0; } /* * subtract out electrostatic contribution from excluded pairs * if pair is 1-4 with scaled1-4 in effect, must subtract out * the scaled difference */ void subtract_excl(Force *force) { MD_Dvec r_ij; MD_Dvec *pos = force->data->pos; MD_Dvec *f = force->data->f; double rlen2; MD_Int **excllist = force->excllist; MD_Int **scaled14 = force->scaled14; MD_Int *excl; MD_Int natoms = force->data->natoms; MD_Int i, j; for (j = 0; j < natoms; j++) { for (excl = excllist[j]; *excl < j; excl++) { i = *excl; r_ij.x = pos[j].x - pos[i].x; r_ij.y = pos[j].y - pos[i].y; r_ij.z = pos[j].z - pos[i].z; rlen2 = r_ij.x * r_ij.x + r_ij.y * r_ij.y + r_ij.z * r_ij.z; subtract_elecpair(force, i, j, rlen2, &r_ij, &(f[i]), &(f[j])); } } if (scaled14 == NULL || force->data->scaling14 == 1.0) return; for (j = 0; j < natoms; j++) { for (excl = scaled14[j]; *excl < j; excl++) { i = *excl; r_ij.x = pos[j].x - pos[i].x; r_ij.y = pos[j].y - pos[i].y; r_ij.z = pos[j].z - pos[i].z; rlen2 = r_ij.x * r_ij.x + r_ij.y * r_ij.y + r_ij.z * r_ij.z; subtract_elec14pair(force, i, j, rlen2, &r_ij, &(f[i]), &(f[j])); } } } /* * evaluate nonbonded forces using box neighbor lists */ void compute_atombox(Force *force, FunptrComputePair compute_pair) { MD_Dvec *pos = force->data->pos; MD_Dvec *f = force->data->f; double *energy = force->data->energy; const double cutoff2 = force->cutoff2; const MD_Int numbox = force->numbox; MD_Dvec r_ij; double rlen2; AtomBox *box = force->atombox; MD_Int i, j, k, ni, nj, nk, nn; /* evaluate nonbonded forces in each box */ for (k = 0; k < numbox; k++) { /* atoms in this box */ for (ni = 1; ni < box[k].num; ni++) { for (nj = 0; nj < ni; nj++) { i = box[k].atomid[ni]; j = box[k].atomid[nj]; r_ij.x = pos[j].x - pos[i].x; r_ij.y = pos[j].y - pos[i].y; r_ij.z = pos[j].z - pos[i].z; rlen2 = r_ij.x * r_ij.x + r_ij.y * r_ij.y + r_ij.z * r_ij.z; if (rlen2 < cutoff2) { compute_pair(force, i, j, rlen2, &r_ij, &(f[i]), &(f[j])); } } } /* atoms between this box and each of its neighbors */ for (nn = 0; nn < box[k].numnbrs; nn++) { nk = box[k].nbr[nn]; for (ni = 0; ni < box[k].num; ni++) { for (nj = 0; nj < box[nk].num; nj++) { i = box[k].atomid[ni]; j = box[nk].atomid[nj]; r_ij.x = pos[j].x - pos[i].x; r_ij.y = pos[j].y - pos[i].y; r_ij.z = pos[j].z - pos[i].z; rlen2 = r_ij.x * r_ij.x + r_ij.y * r_ij.y + r_ij.z * r_ij.z; if (rlen2 < cutoff2) { compute_pair(force, i, j, rlen2, &r_ij, &(f[i]), &(f[j])); } } } } } } /* * evaluate VDW and electrostatic forces between a pair of atoms * (assume that the pair of atoms is within the cutoff distance) */ void compute_nbpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j) { MD_Atom *atom = force->data->atom; const MD_Int natomprms = force->data->natomprms; const double cutoff2 = force->cutoff2; const double inv_cutoff2 = force->inv_cutoff2; const double switchdist2 = force->switchdist2; const double epsilon14 = force->data->scaling14; MD_Int **excllist = force->excllist; MD_Int **scaled14 = force->scaled14; const MD_Int *excl; MD_Dvec f_nb; double rlen, inv_rlen, inv_rlen2, inv_rlen6, inv_rlen12; double c_elec, sw_elec, c_nb; double aterm_vdw, bterm_vdw, efac_vdw, ffac_vdw; double swa, swb, swv, dswv; const double *entry; /* check to see if this pair is excluded */ for (excl = excllist[i]; *excl < j; excl++) ; if (j == *excl) return; /* compute simple functions of rlen2 */ rlen = sqrt(rlen2); inv_rlen = 1.0 / rlen; inv_rlen2 = inv_rlen * inv_rlen; inv_rlen6 = inv_rlen2 * inv_rlen2 * inv_rlen2; inv_rlen12 = inv_rlen6 * inv_rlen6; /* electrostatic computations */ sw_elec = 1.0 - rlen2 * inv_cutoff2; c_elec = force->elec_const * atom[i].q * atom[j].q * inv_rlen * sw_elec; /* choose correct van der Waals parameters */ entry = force->vdwtable + 4 * (atom[i].param * natomprms + atom[j].param); /* determine whether this is a 1-4 interaction */ if (scaled14 == NULL) { aterm_vdw = entry[A] * inv_rlen12; bterm_vdw = entry[B] * inv_rlen6; } else { for (excl = scaled14[i]; *excl < j; excl++) ; if (j != *excl) { aterm_vdw = entry[A] * inv_rlen12; bterm_vdw = entry[B] * inv_rlen6; } else { aterm_vdw = entry[A_14] * inv_rlen12; bterm_vdw = entry[B_14] * inv_rlen6; c_elec *= epsilon14; } } /* van der Waals computations */ efac_vdw = aterm_vdw - bterm_vdw; ffac_vdw = (12.0 * aterm_vdw - 6.0 * bterm_vdw) * inv_rlen2; if (force->is_switching && rlen2 > switchdist2) { swa = cutoff2 - rlen2; swb = swa * (cutoff2 + 2.0 * rlen2 - 3.0 * switchdist2); swv = swa * swb * force->sw_denom; dswv = (swa * swa - swb) * force->sw_denom_times_four; ffac_vdw = ffac_vdw * swv - efac_vdw * dswv; efac_vdw *= swv; } force->pe_vdw += efac_vdw; /* finish electrostatic computations */ force->pe_elec += c_elec * sw_elec; c_elec *= 3.0 * inv_cutoff2 + inv_rlen2; /* find and accumulate nonbonded forces */ c_nb = c_elec + ffac_vdw; f_nb.x = c_nb * r_ij->x; f_nb.y = c_nb * r_ij->y; f_nb.z = c_nb * r_ij->z; f_i->x -= f_nb.x; f_i->y -= f_nb.y; f_i->z -= f_nb.z; f_j->x += f_nb.x; f_j->y += f_nb.y; f_j->z += f_nb.z; } /* * compute VDW forces between pair of atoms * check exclusion list, use switching function */ void compute_vdwpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j) { MD_Atom *atom = force->data->atom; const MD_Int natomprms = force->data->natomprms; const double cutoff2 = force->cutoff2; const double switchdist2 = force->switchdist2; MD_Int **excllist = force->excllist; MD_Int **scaled14 = force->scaled14; const MD_Int *excl; MD_Dvec f_vdw; double inv_rlen2, inv_rlen6, inv_rlen12; double aterm_vdw, bterm_vdw, efac_vdw, ffac_vdw; double swa, swb, swv, dswv; const double *entry; /* check to see if this pair is excluded */ for (excl = excllist[i]; *excl < j; excl++) ; if (j == *excl) return; /* compute simple functions of rlen2 */ inv_rlen2 = 1.0 / rlen2; inv_rlen6 = inv_rlen2 * inv_rlen2 * inv_rlen2; inv_rlen12 = inv_rlen6 * inv_rlen6; /* choose correct van der Waals parameters */ entry = force->vdwtable + 4 * (atom[i].param * natomprms + atom[j].param); /* determine whether this is a 1-4 interaction */ if (scaled14 == NULL) { aterm_vdw = entry[A] * inv_rlen12; bterm_vdw = entry[B] * inv_rlen6; } else { for (excl = scaled14[i]; *excl < j; excl++) ; if (j != *excl) { aterm_vdw = entry[A] * inv_rlen12; bterm_vdw = entry[B] * inv_rlen6; } else { aterm_vdw = entry[A_14] * inv_rlen12; bterm_vdw = entry[B_14] * inv_rlen6; } } /* van der Waals computations */ efac_vdw = aterm_vdw - bterm_vdw; ffac_vdw = (12.0 * aterm_vdw - 6.0 * bterm_vdw) * inv_rlen2; if (force->is_switching && rlen2 > switchdist2) { swa = cutoff2 - rlen2; swb = swa * (cutoff2 + 2.0 * rlen2 - 3.0 * switchdist2); swv = swa * swb * force->sw_denom; dswv = (swa * swa - swb) * force->sw_denom_times_four; ffac_vdw = ffac_vdw * swv - efac_vdw * dswv; efac_vdw *= swv; } force->pe_vdw += efac_vdw; /* find and accumulate nonbonded forces */ f_vdw.x = ffac_vdw * r_ij->x; f_vdw.y = ffac_vdw * r_ij->y; f_vdw.z = ffac_vdw * r_ij->z; f_i->x -= f_vdw.x; f_i->y -= f_vdw.y; f_i->z -= f_vdw.z; f_j->x += f_vdw.x; f_j->y += f_vdw.y; f_j->z += f_vdw.z; } /* * compute electrostatic forces between pair of atoms * DO NOT use shifting function, DO NOT check exclusion list */ void compute_elecpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j) { MD_Atom *atom = force->data->atom; MD_Dvec f_elec; double rlen, inv_rlen; double c_elec; /* compute simple functions of rlen2 */ rlen = sqrt(rlen2); inv_rlen = 1.0 / rlen; /* electrostatic computations */ c_elec = force->elec_const * atom[i].q * atom[j].q * inv_rlen; force->pe_elec += c_elec; c_elec *= inv_rlen * inv_rlen; /* find and accumulate electrostatic forces */ f_elec.x = c_elec * r_ij->x; f_elec.y = c_elec * r_ij->y; f_elec.z = c_elec * r_ij->z; f_i->x -= f_elec.x; f_i->y -= f_elec.y; f_i->z -= f_elec.z; f_j->x += f_elec.x; f_j->y += f_elec.y; f_j->z += f_elec.z; } /* * subtract the electrostatic potential and forces due to this atom pair * DO NOT use shifting function */ void subtract_elecpair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j) { MD_Atom *atom = force->data->atom; MD_Dvec f_elec; double rlen, inv_rlen; double c_elec; /* compute simple functions of rlen2 */ rlen = sqrt(rlen2); inv_rlen = 1.0 / rlen; /* negate electrostatic computations */ c_elec = -force->elec_const * atom[i].q * atom[j].q * inv_rlen; force->pe_elec += c_elec; c_elec *= inv_rlen * inv_rlen; /* find and accumulate opposite electrostatic forces */ f_elec.x = c_elec * r_ij->x; f_elec.y = c_elec * r_ij->y; f_elec.z = c_elec * r_ij->z; f_i->x -= f_elec.x; f_i->y -= f_elec.y; f_i->z -= f_elec.z; f_j->x += f_elec.x; f_j->y += f_elec.y; f_j->z += f_elec.z; } /* * subtract the regular electrostatic potential and forces * due to this atom pair and add in those for a 1-4 interaction * DO NOT use shifting function */ void subtract_elec14pair(Force *force, MD_Int i, MD_Int j, double rlen2, MD_Dvec *r_ij, MD_Dvec *f_i, MD_Dvec *f_j) { MD_Atom *atom = force->data->atom; MD_Dvec f_elec; double rlen, inv_rlen; double c_elec; /* compute simple functions of rlen2 */ rlen = sqrt(rlen2); inv_rlen = 1.0 / rlen; /* negate electrostatic computations */ c_elec = (force->data->scaling14 - 1.0) * force->elec_const * atom[i].q * atom[j].q * inv_rlen; force->pe_elec += c_elec; c_elec *= inv_rlen * inv_rlen; /* find and accumulate opposite electrostatic forces */ f_elec.x = c_elec * r_ij->x; f_elec.y = c_elec * r_ij->y; f_elec.z = c_elec * r_ij->z; f_i->x -= f_elec.x; f_i->y -= f_elec.y; f_i->z -= f_elec.z; f_j->x += f_elec.x; f_j->y += f_elec.y; f_j->z += f_elec.z; }