--- branches/development/src/nonbonded/GB.cpp 2010/10/02 19:53:32 1502 +++ branches/development/src/nonbonded/GB.cpp 2012/03/10 04:21:44 1686 @@ -36,7 +36,8 @@ * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -222,31 +223,45 @@ namespace OpenMD { dw2 = 1.0; } - GBInteractionData mixer; + GBInteractionData mixer1, mixer2; // Cleaver paper uses sqrt of squares to get sigma0 for // mixed interactions. - mixer.sigma0 = sqrt(d1*d1 + d2*d2); - mixer.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2); - mixer.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1); - mixer.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) / + mixer1.sigma0 = sqrt(d1*d1 + d2*d2); + mixer1.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2); + mixer1.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1); + mixer1.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) / ((l2*l2 + d1*d1) * (l1*l1 + d2*d2)); + + mixer2.sigma0 = mixer1.sigma0; + // xa2 and xai2 for j-i pairs are reversed from the same i-j pairing. + // Swapping the particles reverses the anisotropy parameters: + mixer2.xa2 = mixer1.xai2; + mixer2.xai2 = mixer1.xa2; + mixer2.x2 = mixer1.x2; // assumed LB mixing rules for now: - mixer.dw = 0.5 * (dw1 + dw2); - mixer.eps0 = sqrt(e1 * e2); + mixer1.dw = 0.5 * (dw1 + dw2); + mixer1.eps0 = sqrt(e1 * e2); + + mixer2.dw = mixer1.dw; + mixer2.eps0 = mixer1.eps0; RealType er = sqrt(er1 * er2); - RealType ermu = pow(er,(1.0 / mu_)); + RealType ermu = pow(er, (RealType(1.0) / mu_)); RealType xp = (1.0 - ermu) / (1.0 + ermu); RealType ap2 = 1.0 / (1.0 + ermu); - mixer.xp2 = xp * xp; - mixer.xpap2 = xp * ap2; - mixer.xpapi2 = xp / ap2; + mixer1.xp2 = xp * xp; + mixer1.xpap2 = xp * ap2; + mixer1.xpapi2 = xp / ap2; + mixer2.xp2 = mixer1.xp2; + mixer2.xpap2 = mixer1.xpap2; + mixer2.xpapi2 = mixer1.xpapi2; + // only add this pairing if at least one of the atoms is a Gay-Berne atom if (atomType->isGayBerne() || atype2->isGayBerne()) { @@ -255,20 +270,19 @@ namespace OpenMD { key1 = make_pair(atomType, atype2); key2 = make_pair(atype2, atomType); - MixingMap[key1] = mixer; + MixingMap[key1] = mixer1; if (key2 != key1) { - MixingMap[key2] = mixer; + MixingMap[key2] = mixer2; } } } } - void GB::calcForce(InteractionData idat) { + void GB::calcForce(InteractionData &idat) { if (!initialized_) initialize(); - pair key = make_pair(idat.atype1, idat.atype2); - GBInteractionData mixer = MixingMap[key]; + GBInteractionData mixer = MixingMap[idat.atypes]; RealType sigma0 = mixer.sigma0; RealType dw = mixer.dw; @@ -280,26 +294,26 @@ namespace OpenMD { RealType xpap2 = mixer.xpap2; RealType xpapi2 = mixer.xpapi2; - Vector3d ul1 = idat.A1.getRow(2); - Vector3d ul2 = idat.A2.getRow(2); + Vector3d ul1 = idat.A1->getRow(2); + Vector3d ul2 = idat.A2->getRow(2); RealType a, b, g; - bool i_is_LJ = idat.atype1->isLennardJones(); - bool j_is_LJ = idat.atype2->isLennardJones(); - + bool i_is_LJ = idat.atypes.first->isLennardJones(); + bool j_is_LJ = idat.atypes.second->isLennardJones(); + if (i_is_LJ) { a = 0.0; ul1 = V3Zero; } else { - a = dot(idat.d, ul1); + a = dot(*(idat.d), ul1); } if (j_is_LJ) { b = 0.0; ul2 = V3Zero; } else { - b = dot(idat.d, ul2); + b = dot(*(idat.d), ul2); } if (i_is_LJ || j_is_LJ) @@ -307,8 +321,8 @@ namespace OpenMD { else g = dot(ul1, ul2); - RealType au = a / idat.rij; - RealType bu = b / idat.rij; + RealType au = a / *(idat.rij); + RealType bu = b / *(idat.rij); RealType au2 = au * au; RealType bu2 = bu * bu; @@ -321,7 +335,7 @@ namespace OpenMD { RealType e1 = 1.0 / sqrt(1.0 - x2*g2); RealType e2 = 1.0 - Hp; RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_); - RealType BigR = dw*sigma0 / (idat.rij - sigma + dw*sigma0); + RealType BigR = dw*sigma0 / (*(idat.rij) - sigma + dw*sigma0); RealType R3 = BigR*BigR*BigR; RealType R6 = R3*R3; @@ -329,16 +343,19 @@ namespace OpenMD { RealType R12 = R6*R6; RealType R13 = R6*R7; - RealType U = idat.vdwMult * 4.0 * eps * (R12 - R6); + RealType U = *(idat.vdwMult) * 4.0 * eps * (R12 - R6); RealType s3 = sigma*sigma*sigma; RealType s03 = sigma0*sigma0*sigma0; - RealType pref1 = - idat.vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * idat.rij); + RealType pref1 = - *(idat.vdwMult) * 8.0 * eps * mu_ * (R12 - R6) / + (e2 * *(idat.rij)); - RealType pref2 = idat.vdwMult * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /(dw*idat.rij*s03); + RealType pref2 = *(idat.vdwMult) * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) / + (dw* *(idat.rij) * s03); - RealType dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*idat.rij/s3 + H)); + RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 * + *(idat.rij) / s3 + H)); RealType dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0 - xp2 * g2) + pref2 * (xa2 * au - x2 *bu*g) / (1.0 - x2 * g2); @@ -350,21 +367,47 @@ namespace OpenMD { + 8.0 * eps * mu_ * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) / (1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) * (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03); - - Vector3d rhat = idat.d / idat.rij; - Vector3d rxu1 = cross(idat.d, ul1); - Vector3d rxu2 = cross(idat.d, ul2); + Vector3d rhat = *(idat.d) / *(idat.rij); + Vector3d rxu1 = cross(*(idat.d), ul1); + Vector3d rxu2 = cross(*(idat.d), ul2); Vector3d uxu = cross(ul1, ul2); - - idat.pot += U*idat.sw; - idat.f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2; - idat.t1 += dUda * rxu1 - dUdg * uxu; - idat.t2 += dUdb * rxu2 - dUdg * uxu; - idat.vpair += U*idat.sw; + (*(idat.pot))[VANDERWAALS_FAMILY] += U * *(idat.sw); + *(idat.f1) += (dUdr * rhat + dUda * ul1 + dUdb * ul2) * *(idat.sw); + *(idat.t1) += (dUda * rxu1 - dUdg * uxu) * *(idat.sw); + *(idat.t2) += (dUdb * rxu2 + dUdg * uxu) * *(idat.sw); + *(idat.vpair) += U; + return; } + + RealType GB::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); + + RealType cut = 0.0; + + if (atypes.first->isGayBerne()) { + GayBerneParam gb1 = getGayBerneParam(atypes.first); + RealType d1 = gb1.GB_d; + RealType l1 = gb1.GB_l; + // sigma is actually sqrt(2)*l for prolate ellipsoids + cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1)); + } else if (atypes.first->isLennardJones()) { + cut = max(cut, RealType(2.5) * getLJSigma(atypes.first)); + } + + if (atypes.second->isGayBerne()) { + GayBerneParam gb2 = getGayBerneParam(atypes.second); + RealType d2 = gb2.GB_d; + RealType l2 = gb2.GB_l; + cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2)); + } else if (atypes.second->isLennardJones()) { + cut = max(cut, RealType(2.5) * getLJSigma(atypes.second)); + } + + return cut; + } }