The Art of Molecular Dynamics Simulation V2 源码

/* [[pr_13_3 - fast-multipole method]] */ /********************************************************************* (C) 2004 D. C. Rapaport This software is copyright material accompanying the book "The Art of Molecular Dynamics Simulation", 2nd edition, by D. C. Rapaport, published by Cambridge University Press (2004). **********************************************************************/ #include "in_mddefs.h" typedef struct { VecR r, rv, ra; real chg; } Mol; Mol *mol; VecR region, vSum; VecI initUcell; real deltaT, density, rCut, temperature, timeNow, uSum, velMag, vvSum; Prop kinEnergy, potEnergy, totEnergy; int moreCycles, nMol, stepAvg, stepCount, stepEquil, stepLimit; VecI cells; int *cellList; real dispHi, rNebrShell; int *nebrTab, nebrNow, nebrTabFac, nebrTabLen, nebrTabMax; real **histRdf, **cumRdf, rangeRdf; int countRdf, limitRdf, sizeHistRdf, stepRdf; real kinEnInitSum; int stepInitlzTemp; MpCell **mpCell; VecR cellWid; VecI mpCells; real chargeMag; int *mpCellList, curCellsEdge, curLevel, maxCellsEdge, maxLevel, maxOrd, wellSep; NameList nameList[] = { NameR (chargeMag), NameR (deltaT), NameR (density), NameI (initUcell), NameI (limitRdf), NameI (maxLevel), NameI (nebrTabFac), NameR (rangeRdf), NameR (rNebrShell), NameI (sizeHistRdf), NameI (stepAvg), NameI (stepEquil), NameI (stepInitlzTemp), NameI (stepLimit), NameI (stepRdf), NameR (temperature), NameI (wellSep), }; int main (int argc, char **argv) { GetNameList (argc, argv); PrintNameList (stdout); SetParams (); SetupJob (); moreCycles = 1; while (moreCycles) { SingleStep (); if (stepCount >= stepLimit) moreCycles = 0; } } void SingleStep () { ++ stepCount; timeNow = stepCount * deltaT; LeapfrogStep (1); if (nebrNow) { nebrNow = 0; dispHi = 0.; BuildNebrList (); } ComputeForces (); MultipoleCalc (); ComputeWallForces (); ApplyThermostat (); LeapfrogStep (2); EvalProps (); if (stepCount < stepEquil) AdjustInitTemp (); AccumProps (1); if (stepCount % stepAvg == 0) { AccumProps (2); PrintSummary (stdout); AccumProps (0); } if (stepCount >= stepEquil && (stepCount - stepEquil) % stepRdf == 0) EvalRdf (); } void SetupJob () { AllocArrays (); stepCount = 0; InitCoords (); InitVels (); InitAccels (); InitCharges (); AccumProps (0); nebrNow = 1; kinEnInitSum = 0.; } void SetParams () { rCut = pow (2., 1./6.); VSCopy (region, 1. / pow (density, 1./3.), initUcell); nMol = VProd (initUcell); velMag = sqrt (NDIM * (1. - 1. / nMol) * temperature); VSCopy (cells, 1. / (rCut + rNebrShell), region); nebrTabMax = nebrTabFac * nMol; maxOrd = MAX_MPEX_ORD; } void AllocArrays () { int k, n; AllocMem (mol, nMol, Mol); AllocMem (cellList, VProd (cells) + nMol, int); AllocMem (nebrTab, 2 * nebrTabMax, int); AllocMem (mpCell, maxLevel + 1, MpCell *); maxCellsEdge = 2; for (n = 2; n <= maxLevel; n ++) { maxCellsEdge *= 2; VSetAll (mpCells, maxCellsEdge); AllocMem (mpCell[n], VProd (mpCells), MpCell); } AllocMem (mpCellList, nMol + VProd (mpCells), int); AllocMem2 (histRdf, 2, sizeHistRdf, real); AllocMem2 (cumRdf, 2, sizeHistRdf, real); } void BuildNebrList () { VecR dr, invWid, rs; VecI cc, m1v, m2v, vOff[] = OFFSET_VALS; real rrNebr; int c, j1, j2, m1, m1x, m1y, m1z, m2, n, offset; rrNebr = Sqr (rCut + rNebrShell); VDiv (invWid, cells, region); for (n = nMol; n < nMol + VProd (cells); n ++) cellList[n] = -1; DO_MOL { VSAdd (rs, mol[n].r, 0.5, region); VMul (cc, rs, invWid); c = VLinear (cc, cells) + nMol; cellList[n] = cellList[c]; cellList[c] = n; } nebrTabLen = 0; for (m1z = 0; m1z < cells.z; m1z ++) { for (m1y = 0; m1y < cells.y; m1y ++) { for (m1x = 0; m1x < cells.x; m1x ++) { VSet (m1v, m1x, m1y, m1z); m1 = VLinear (m1v, cells) + nMol; for (offset = 0; offset < N_OFFSET; offset ++) { VAdd (m2v, m1v, vOff[offset]); if (m2v.x < 0 || m2v.x >= cells.x || m2v.y < 0 || m2v.y >= cells.y || m2v.z >= cells.z) continue; m2 = VLinear (m2v, cells) + nMol; DO_CELL (j1, m1) { DO_CELL (j2, m2) { if (m1 != m2 || j2 < j1) { VSub (dr, mol[j1].r, mol[j2].r); if (VLenSq (dr) < rrNebr) { if (nebrTabLen >= nebrTabMax) ErrExit (ERR_TOO_MANY_NEBRS); nebrTab[2 * nebrTabLen] = j1; nebrTab[2 * nebrTabLen + 1] = j2; ++ nebrTabLen; } } } } } } } } } void ComputeForces () { VecR dr; real fcVal, rr, rrCut, rri, rri3, uVal; int j1, j2, n; rrCut = Sqr (rCut); DO_MOL VZero (mol[n].ra); uSum = 0.; for (n = 0; n < nebrTabLen; n ++) { j1 = nebrTab[2 * n]; j2 = nebrTab[2 * n + 1]; VSub (dr, mol[j1].r, mol[j2].r); rr = VLenSq (dr); if (rr < rrCut) { rri = 1. / rr; rri3 = Cube (rri); fcVal = 48. * rri3 * (rri3 - 0.5) * rri; uVal = 4. * rri3 * (rri3 - 1.) + 1.; VVSAdd (mol[j1].ra, fcVal, dr); VVSAdd (mol[j2].ra, - fcVal, dr); uSum += uVal; } } } void MultipoleCalc () { int j, k, m1; VSetAll (mpCells, maxCellsEdge); AssignMpCells (); VDiv (cellWid, region, mpCells); EvalMpCell (); curCellsEdge = maxCellsEdge; for (curLevel = maxLevel - 1; curLevel >= 2; curLevel --) { curCellsEdge /= 2; VSetAll (mpCells, curCellsEdge); VDiv (cellWid, region, mpCells); CombineMpCell (); } for (m1 = 0; m1 < 64; m1 ++) { for (j = 0; j <= maxOrd; j ++) { for (k = 0; k <= j; k ++) { mpCell[2][m1].me.c(j, k) = 0.; mpCell[2][m1].me.s(j, k) = 0.; } } } curCellsEdge = 2; for (curLevel = 2; curLevel <= maxLevel; curLevel ++) { curCellsEdge *= 2; VSetAll (mpCells, curCellsEdge); VDiv (cellWid, region, mpCells); GatherWellSepLo (); if (curLevel < maxLevel) PropagateCellLo (); } ComputeFarCellInt (); ComputeNearCellInt (); } void AssignMpCells () { VecR invWid, rs; VecI cc; int c, n; VDiv (invWid, mpCells, region); for (n = nMol; n < nMol + VProd (mpCells); n ++) mpCellList[n] = -1; DO_MOL { VSAdd (rs, mol[n].r, 0.5, region); VMul (cc, rs, invWid); c = VLinear (cc, mpCells) + nMol; mpCellList[n] = mpCellList[c]; mpCellList[c] = n; } } #define DO_MP_CELL(j, m) \ for (j = mpCellList[m + nMol]; j >= 0; j = mpCellList[j]) void EvalMpCell () { MpTerms le; VecR cMid, dr; VecI m1v; int j, j1, k, m1, m1x, m1y, m1z; for (m1z = 0; m1z < mpCells.z; m1z ++) { for (m1y = 0; m1y < mpCells.y; m1y ++) { for (m1x = 0; m1x < mpCells.x; m1x ++) { VSet (m1v, m1x, m1y, m1z); m1 = VLinear (m1v, mpCells); mpCell[maxLevel][m1].occ = 0; for (j = 0; j <= maxOrd; j ++) { for (k = 0; k <= j; k ++) { mpCell[maxLevel][m1].le.c(j, k) = 0.; mpCell[maxLevel][m1].le.s(j, k) = 0.; } } if (mpCellList[m1 + nMol] >= 0) { VAddCon (cMid, m1v, 0.5); VMul (cMid, cMid, cellWid); VVSAdd (cMid, - 0.5, region); DO_MP_CELL (j1, m1) { ++ mpCell[maxLevel][m1].occ; VSub (dr, mol[j1].r, cMid); EvalMpL (&le, &dr, maxOrd); for (j = 0; j <= maxOrd; j ++) { for (k = 0; k <= j; k ++) { mpCell[maxLevel][m1].le.c(j, k) += mol[j1].chg * le.c(j, k); mpCell[maxLevel][m1].le.s(j, k) += mol[j1].chg * le.s(j, k); } } } } } } } } void CombineMpCell () { MpTerms le, le2; VecR rShift; VecI m1v, m2v, mpCellsN; int iDir, j, k, m1, m1x, m1y, m1z, m2; VSCopy (mpCellsN, 2, mpCells);