/* * force.c * * High level force routines, build exclusion lists and vdw param matrix. */ #include #include #include #include "mdutil.h" #include "force.h" #include "fbond.h" #include "fnonbond.h" /* * prototypes for internal functions */ static double compute_sphericalbc(Force *force); static MD_Errcode build_vdwtable(Force *force); static MD_Errcode build_exclusions(Force *force); static void build_exclusions_freemem(MD_Int natoms, MD_Int **exclx, MD_Int *lenx, MD_Int **excl12, MD_Int *len12, MD_Int **excl13, MD_Int *len13, MD_Int **excl14, MD_Int **scaled14, MD_Int *accum, MD_Int *dest); static void sort(MD_Int *list, MD_Int len); static MD_Int merge(MD_Int *dest, const MD_Int *src1, const MD_Int *src2, MD_Int n); /* * initialize the force object * * assumes that force is already allocated and memory cleared * data must already be initialized * * allocate and initialize nonbonded exclusion lists * * allocate and initialize van der Waals parameters table */ MD_Errcode deven_force_init(Force *force, Data *data) { /* set pointer back to engine */ force->data = data; /* allocate and initialize nonbonded force structures */ if (build_vdwtable(force)) return MD_FAIL; if (build_exclusions(force)) return MD_FAIL; /* calculate nonbonded force constants */ force->elec_const = COULOMB / data->dielectric; force->cutoff2 = data->cutoff * data->cutoff; force->inv_cutoff2 = 1.0 / force->cutoff2; force->switchdist2 = data->switchdist * data->switchdist; force->sw_denom = 1.0 / (force->cutoff2 - force->switchdist2); force->sw_denom = force->sw_denom * force->sw_denom * force->sw_denom; force->sw_denom_times_four = 4.0 * force->sw_denom; force->is_switching = data->is_switching; /* init parameters for spherical boundary conditions */ /* make sure that "r1" is radius of smaller sphere, swap if not */ force->is_sphericalbc = data->is_sphericalbc; force->center = data->sphericalbc_center; if (data->sphericalbc_k2 == 0.0 || data->sphericalbc_r1 < data->sphericalbc_r2) { force->r1 = data->sphericalbc_r1; force->k1 = data->sphericalbc_k1 * INV_KCAL_PER_MOL; force->p1 = data->sphericalbc_exp1; force->r2 = data->sphericalbc_r2; force->k2 = data->sphericalbc_k2 * INV_KCAL_PER_MOL; force->p2 = data->sphericalbc_exp2; } else { force->r1 = data->sphericalbc_r2; force->k1 = data->sphericalbc_k2 * INV_KCAL_PER_MOL; force->p1 = data->sphericalbc_exp2; force->r2 = data->sphericalbc_r1; force->k2 = data->sphericalbc_k1 * INV_KCAL_PER_MOL; force->p2 = data->sphericalbc_exp1; } if (deven_nonbonded_init(force)) return MD_FAIL; return 0; } /* * destroy the force object * * free all memory allocated and clear memory used by force object */ void deven_force_destroy(Force *force) { MD_Int natoms = 0; MD_Int i; if (force->data) natoms = force->data->natoms; if (force->excllist) { for (i = 0; i < natoms; i++) free(force->excllist[i]); } if (force->scaled14) { for (i = 0; i < natoms; i++) free(force->scaled14[i]); } free(force->excllist); free(force->scaled14); free(force->vdwtable); deven_nonbonded_destroy(force); memset(force, 0, sizeof(Force)); } /* * compute the force and the potential */ MD_Errcode deven_force_compute(Force *force) { MD_Double *energy = force->data->energy; MD_Double pe = 0.0; /* zero the force vector before accumulating */ memset(force->data->f, 0, force->data->natoms * sizeof(MD_Dvec)); /* compute bonded forces */ pe += energy[PE_BOND] = deven_compute_bonds(force); pe += energy[PE_ANGLE] = deven_compute_angles(force); pe += energy[PE_DIHED] = deven_compute_dihedrals(force); pe += energy[PE_IMPR] = deven_compute_impropers(force); /* compute nonbonded forces */ force->pe_elec = force->pe_vdw = 0.0; if (force->compute_nonbonded(force)) return MD_FAIL; pe += (energy[PE_ELEC]=force->pe_elec) + (energy[PE_VDW]=force->pe_vdw); /* spherical boundary conditions */ if (force->is_sphericalbc == ON) { pe += energy[PE_BC] = compute_sphericalbc(force); } /* store accumulated potential energy */ energy[PE] = pe; return 0; } /* * compute spherical boundary conditions */ double compute_sphericalbc(Force *force) { MD_Dvec *pos = force->data->pos; MD_Dvec *f = force->data->f; const MD_Int natoms = force->data->natoms; const MD_Dvec center = force->center; const double r1 = force->r1; const double k1 = force->k1; const MD_Int p1 = force->p1; const double r2 = force->r2; const double k2 = force->k2; const MD_Int p2 = force->p2; const MD_Int r1_sq = r1 * r1; const MD_Int r2_sq = r2 * r2; double dx, dy, dz, dist, dist2, d1, d2, fac1, fac2; double pe_sphbc = 0.0, f_sphbc; MD_Int i, j; for (i = 0; i < natoms; i++) { dx = center.x - pos[i].x; dy = center.y - pos[i].y; dz = center.z - pos[i].z; dist2 = dx * dx + dy * dy + dz * dz; if (dist2 <= r1_sq) continue; dist = sqrt(dist2); d1 = dist - r1; fac1 = k1; for (j = 1; j < p1; j++) fac1 *= d1; if (k2 == 0.0 || dist2 <= r2_sq) { pe_sphbc += fac1 * d1; f_sphbc = p1 * fac1 / dist; } else { d2 = dist - r2; fac2 = k2; for (j = 1; j < p2; j++) fac2 *= d2; pe_sphbc += fac1 * d1 + fac2 * d2; f_sphbc = (p1 * fac1 + p2 * fac2) / dist; } f[i].x += f_sphbc * dx; f[i].y += f_sphbc * dy; f[i].z += f_sphbc * dz; } return pe_sphbc; } /* * build the van der Waals parameter table * * table is a "square" symmetric matrix, dimension (natomprms * natomprms) * each entry of matrix contains A, B, A_14, B_14 parameters * matrix is indexed by atom "types" (0..natomprms-1) * * vdwtable is stored as one-dimensional array * index for (i,j) atom pair interaction is: 4 * (i * natomprms + j) */ MD_Errcode build_vdwtable(Force *force) { MD_Atom_Param *atomprm = force->data->atomprm; const MD_Int natomprms = force->data->natomprms; MD_Nbfix_Param *nbfixprm = force->data->nbfixprm; const MD_Int nnbfixprms = force->data->nnbfixprms; double *vdwtable, *ij_entry, *ji_entry; double neg_emin, rmin, neg_emin14, rmin14; MD_Int i, j, k; vdwtable = (double *) malloc(4 * natomprms * natomprms * sizeof(double)); if (vdwtable == NULL) { MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_vdwtable call to malloc"); return MD_FAIL; } /* compute each table entry given separate i and j atom params */ for (i = 0; i < natomprms; i++) { for (j = i; j < natomprms; j++) { ij_entry = vdwtable + 4 * (i * natomprms + j); ji_entry = vdwtable + 4 * (j * natomprms + i); /* compute vdw A and B coefficients for atom type ij interaction */ neg_emin = sqrt(atomprm[i].emin * atomprm[j].emin) * INV_KCAL_PER_MOL; rmin = 0.5 * (atomprm[i].rmin + atomprm[j].rmin); neg_emin14 = sqrt(atomprm[i].emin14 * atomprm[j].emin14) * INV_KCAL_PER_MOL; rmin14 = 0.5 * (atomprm[i].rmin14 + atomprm[j].rmin14); /* raise rmin and rmin14 to 6th power */ rmin *= rmin * rmin; rmin *= rmin; rmin14 *= rmin14 * rmin14; rmin14 *= rmin14; /* set ij entry and its transpose */ ij_entry[A] = ji_entry[A] = neg_emin * rmin * rmin; ij_entry[B] = ji_entry[B] = 2.0 * neg_emin * rmin; ij_entry[A_14] = ji_entry[A_14] = neg_emin14 * rmin14 * rmin14; ij_entry[B_14] = ji_entry[B_14] = 2.0 * neg_emin14 * rmin14; } } /* now go back and update entries for nbfix params */ for (k = 0; k < nnbfixprms; k++) { i = nbfixprm[k].param[0]; j = nbfixprm[k].param[1]; ij_entry = vdwtable + 4 * (i * natomprms + j); ji_entry = vdwtable + 4 * (j * natomprms + i); /* compute vdw A and B coefficients for this fixed type interaction */ neg_emin = -nbfixprm[k].emin * INV_KCAL_PER_MOL; rmin = nbfixprm[k].rmin; neg_emin14 = -nbfixprm[k].emin14 * INV_KCAL_PER_MOL; rmin14 = nbfixprm[k].rmin14; /* raise rmin and rmin14 to 6th power */ rmin *= rmin * rmin; rmin *= rmin; rmin14 *= rmin14 * rmin14; rmin14 *= rmin14; /* set ij entry and its transpose */ ij_entry[A] = ji_entry[A] = neg_emin * rmin * rmin; ij_entry[B] = ji_entry[B] = 2.0 * neg_emin * rmin; ij_entry[A_14] = ji_entry[A_14] = neg_emin14 * rmin14 * rmin14; ij_entry[B_14] = ji_entry[B_14] = 2.0 * neg_emin14 * rmin14; } force->vdwtable = vdwtable; return 0; } /* * build the exclusion lists * * lists are built from MD_Excl and MD_Bond * * algorithm (using set notation): * exclx[i] = { j : there is an explicit exclusion (i,j) } * excl12[i] = { j : there is a bond (i,j) } * excl13[i] = (excl12[i] U ( U_{ j \in excl12[i] } excl12[j] )) \ {i} * excl14[i] = (excl13[i] U ( U_{ j \in excl13[i] } excl12[j] )) \ {i} * scaled14[i] = excl14[i] \ excl13[i] * * force->excllist[i] = exclx[i], if excl_policy == EXCL_NONE * = exclx[i] U excl12[i], if excl_policy == EXCL_12 * = exclx[i] U excl13[i], if excl_policy == EXCL_13 * = exclx[i] U excl14[i], if excl_policy == EXCL_14 * * force->excllist[i] = exclx[i] U excl13[i] * AND * force->scaled14[i] = scaled14[i], if excl_policy == SCALED_14 * * allocate little extra memory * implement by merging sorted exclusion lists * each atom's exclusion array is terminated by INT_MAX sentinel */ MD_Errcode build_exclusions(Force *force) { MD_Int **exclx = NULL, *lenx = NULL; MD_Int **excl12 = NULL, *len12 = NULL; MD_Int **excl13 = NULL, *len13 = NULL; MD_Int **excl14 = NULL; MD_Int **scaled14 = NULL; MD_Int *accum = NULL, *dest = NULL; MD_Int *list; MD_Int len, tmp, i, j, k, ii, jj, kk, atom1, atom2; MD_Int size; /* allocated length of accum and dest arrays */ MD_Int maxsize; /* largest length needed to be held by accum or dest */ MD_Int accumlen; /* used length of accum (accumlen <= maxsize <= size) */ MD_Int natoms = force->data->natoms; MD_Int nexcls = force->data->nexcls; MD_Int nbonds = force->data->nbonds; MD_Excl *excl = force->data->excl; MD_Bond *bond = force->data->bond; MD_Int excl_policy = force->data->excl_policy; /* initialize */ force->excllist = NULL; force->scaled14 = NULL; /* we're done if there are no atoms */ if (natoms == 0) return 0; /* allocate memory for explicit exclusions list */ if ((exclx = (MD_Int **) calloc(natoms, sizeof(MD_Int *))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } if ((lenx = (MD_Int *) calloc(natoms, sizeof(MD_Int))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } /* count number of explicit exclusions for each atom */ for (i = 0; i < nexcls; i++) { if (excl[i].atom[0] != excl[i].atom[1]) { lenx[ excl[i].atom[0] ]++; lenx[ excl[i].atom[1] ]++; } } /* allocate memory for each row of exclx, leave space for sentinel */ for (i = 0; i < natoms; i++) { exclx[i] = (MD_Int *) malloc((lenx[i] + 1) * sizeof(MD_Int)); if (exclx[i] == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } lenx[i] = 0; /* zero this to be length counter */ } /* loop over explicit exclusions to fill in the rows of exclx */ for (i = 0; i < nexcls; i++) { atom1 = excl[i].atom[0]; atom2 = excl[i].atom[1]; exclx[atom1][ lenx[atom1]++ ] = atom2; exclx[atom2][ lenx[atom2]++ ] = atom1; } /* place sentinel at end of each row */ for (i = 0; i < natoms; i++) { exclx[i][ lenx[i] ] = INT_MAX; } /* sort each exclx row */ for (i = 0; i < natoms; i++) { sort(exclx[i], lenx[i]); } /* if we're doing no bond exclusions, we're done */ if (excl_policy == EXCL_NONE) { force->excllist = exclx; build_exclusions_freemem(natoms, NULL, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); return 0; } /* allocate memory for 1-2 exclusions list */ if ((excl12 = (MD_Int **) calloc(natoms, sizeof(MD_Int *))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } if ((len12 = (MD_Int *) calloc(natoms, sizeof(MD_Int))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } /* find the length of each row of excl12 */ for (i = 0; i < nbonds; i++) { len12[ bond[i].atom[0] ]++; len12[ bond[i].atom[1] ]++; } /* allocate memory for each row of excl12 */ /* leave space for explicit exclusion list and sentinel */ /* also determine maxsize */ maxsize = 0; for (i = 0; i < natoms; i++) { if (maxsize < len12[i] + lenx[i] + 1) { maxsize = len12[i] + lenx[i] + 1; } excl12[i] = (MD_Int *) malloc((len12[i] + lenx[i] + 1) * sizeof(MD_Int)); if (excl12[i] == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } len12[i] = 0; /* zero this to be length counter */ } /* loop over bonds to fill in the rows of excl12 */ for (i = 0; i < nbonds; i++) { atom1 = bond[i].atom[0]; atom2 = bond[i].atom[1]; excl12[atom1][ len12[atom1]++ ] = atom2; excl12[atom2][ len12[atom2]++ ] = atom1; } /* place sentinel at end of each row */ for (i = 0; i < natoms; i++) { excl12[i][ len12[i] ] = INT_MAX; } /* initialize accum and dest arrays for merge and swap */ size = 10; while (size < maxsize) size *= 2; if ((accum = (MD_Int *) malloc(size * sizeof(MD_Int))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } if ((dest = (MD_Int *) malloc(size * sizeof(MD_Int))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } /* sort each excl12 row */ for (i = 0; i < natoms; i++) { sort(excl12[i], len12[i]); } /* if we're excluding only 1-2 interactions, we're done */ if (excl_policy == EXCL_12) { /* merge each excl12 row with exclx row */ for (i = 0; i < natoms; i++) { len = merge(dest, exclx[i], excl12[i], i); memcpy(excl12[i], dest, (len + 1) * sizeof(MD_Int)); } force->excllist = excl12; build_exclusions_freemem(natoms, exclx, lenx, NULL, len12, excl13, len13, excl14, scaled14, accum, dest); return 0; } /* allocate memory for 1-3 exclusions list */ if ((excl13 = (MD_Int **) calloc(natoms, sizeof(MD_Int *))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } if ((len13 = (MD_Int *) malloc(natoms * sizeof(MD_Int))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } /* merge the excl12 lists into excl13 lists */ for (i = 0; i < natoms; i++) { memcpy(accum, excl12[i], (len12[i] + 1) * sizeof(MD_Int)); accumlen = len12[i]; for (j = 0; excl12[i][j] < INT_MAX; j++) { k = excl12[i][j]; if (k == i) continue; maxsize = accumlen + len12[k]; if (size <= maxsize + lenx[i]) { do { size *= 2; } while (size <= maxsize + lenx[i]); list = (MD_Int *) realloc(accum, size * sizeof(MD_Int)); if (list == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to realloc"); return MD_FAIL; } accum = list; list = (MD_Int *) realloc(dest, size * sizeof(MD_Int)); if (list == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to realloc"); return MD_FAIL; } dest = list; } accumlen = merge(dest, accum, excl12[k], i); list = accum; accum = dest; dest = list; } excl13[i] = (MD_Int *) malloc((accumlen + lenx[i] + 1) * sizeof(MD_Int)); if (excl13[i] == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } memcpy(excl13[i], accum, (accumlen + 1) * sizeof(MD_Int)); len13[i] = accumlen; } /* if we're excluding 1-2 and 1-3 interactions, we're done */ if (excl_policy == EXCL_13) { /* merge each excl13 row with exclx row */ for (i = 0; i < natoms; i++) { len = merge(dest, exclx[i], excl13[i], i); memcpy(excl13[i], dest, (len + 1) * sizeof(MD_Int)); } force->excllist = excl13; build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, NULL, len13, excl14, scaled14, accum, dest); return 0; } /* allocate memory for 1-4 exclusions list */ if ((excl14 = (MD_Int **) calloc(natoms, sizeof(MD_Int *))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } /* merge the excl13 lists into excl14 lists */ for (i = 0; i < natoms; i++) { memcpy(accum, excl13[i], (len13[i] + 1) * sizeof(MD_Int)); accumlen = len13[i]; for (j = 0; excl13[i][j] < INT_MAX; j++) { k = excl13[i][j]; if (k == i) continue; maxsize = accumlen + len12[k]; if (size <= maxsize + lenx[i]) { do { size *= 2; } while (size <= maxsize + lenx[i]); list = (MD_Int *) realloc(accum, size * sizeof(MD_Int)); if (list == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to realloc"); return MD_FAIL; } accum = list; list = (MD_Int *) realloc(dest, size * sizeof(MD_Int)); if (list == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to realloc"); return MD_FAIL; } dest = list; } accumlen = merge(dest, accum, excl12[k], i); list = accum; accum = dest; dest = list; } excl14[i] = (MD_Int *) malloc((accumlen + lenx[i] + 1) * sizeof(MD_Int)); if (excl14[i] == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } memcpy(excl14[i], accum, (accumlen + 1) * sizeof(MD_Int)); } /* if we're excluding 1-2, 1-3, and 1-4 interactions, we're done */ if (excl_policy == EXCL_14) { /* merge each excl14 row with exclx row */ for (i = 0; i < natoms; i++) { len = merge(dest, exclx[i], excl14[i], i); memcpy(excl14[i], dest, (len + 1) * sizeof(MD_Int)); } force->excllist = excl14; build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, NULL, scaled14, accum, dest); return 0; } /* allocate memory for scaled 1-4 list */ if ((scaled14 = (MD_Int **) calloc(natoms, sizeof(MD_Int *))) == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to calloc"); return MD_FAIL; } /* scaled14 list includes everything in excl14 that is not in excl13 */ for (i = 0; i < natoms; i++) { ii = jj = kk = 0; while (excl14[i][ii] < INT_MAX) { if (excl14[i][ii] != excl13[i][jj]) dest[kk++] = excl14[i][ii++]; else { ii++; jj++; } } dest[kk] = INT_MAX; scaled14[i] = (MD_Int *) malloc((kk + 1) * sizeof(MD_Int)); if (scaled14[i] == NULL) { build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, excl13, len13, excl14, scaled14, accum, dest); MD_ERR(force->data->sim, MD_engine_name(force->data->sim), MD_ERR_MEMALLOC, "build_exclusions call to malloc"); return MD_FAIL; } memcpy(scaled14[i], dest, (kk + 1) * sizeof(MD_Int)); } /* set pointers to the lists we need to keep */ /* merge each excl13 row with exclx row */ for (i = 0; i < natoms; i++) { len = merge(dest, exclx[i], excl13[i], i); memcpy(excl13[i], dest, (len + 1) * sizeof(MD_Int)); } force->excllist = excl13; force->scaled14 = scaled14; /* free the memory that is no longer needed */ build_exclusions_freemem(natoms, exclx, lenx, excl12, len12, NULL, len13, excl14, NULL, accum, dest); return 0; } /* * reclaim extra memory allocated in build_exclusions */ void build_exclusions_freemem(MD_Int natoms, MD_Int **exclx, MD_Int *lenx, MD_Int **excl12, MD_Int *len12, MD_Int **excl13, MD_Int *len13, MD_Int **excl14, MD_Int **scaled14, MD_Int *accum, MD_Int *dest) { MD_Int k; if (exclx != NULL) { for (k = 0; k < natoms; k++) free(exclx[k]); free(exclx); } if (excl12 != NULL) { for (k = 0; k < natoms; k++) free(excl12[k]); free(excl12); } if (excl13 != NULL) { for (k = 0; k < natoms; k++) free(excl13[k]); free(excl13); } if (excl14 != NULL) { for (k = 0; k < natoms; k++) free(excl14[k]); free(excl14); } if (scaled14 != NULL) { for (k = 0; k < natoms; k++) free(scaled14[k]); free(scaled14); } free(lenx); free(len12); free(len13); free(accum); free(dest); } /* * sort an array of integers * (use insertion sort: optimal for short arrays) * * assume INT_MAX sentinel is at end of array */ void sort(MD_Int *list, MD_Int len) { MD_Int i, j, tmp; for (i = len - 2; i >= 0; i--) { tmp = list[i]; j = i; while (tmp > list[j+1]) { list[j] = list[j+1]; j++; } list[j] = tmp; } } /* * merge two sorted source arrays into a destination array, * keeping destination sorted and deleting duplicate entries * and excluding n from being merged (used for the self entry) * * assume destination array has enough space * assume each source array is terminated by sentinel INT_MAX * add terminating sentinel INT_MAX to destination array * * return length of destination (not including sentinel) */ MD_Int merge(MD_Int *dest, const MD_Int *src1, const MD_Int *src2, MD_Int n) { MD_Int i = 0, j = 0, k = 0; while (src1[i] < INT_MAX || src2[j] < INT_MAX) { if (src1[i] == n) i++; else if (src2[j] == n) j++; else if (src1[i] < src2[j]) dest[k++] = src1[i++]; else if (src1[i] > src2[j]) dest[k++] = src2[j++]; else { dest[k++] = src1[i++]; j++; } } dest[k] = INT_MAX; return k; }