--- branches/development/src/nonbonded/Sticky.cpp 2010/07/28 19:52:00 1485 +++ branches/development/src/nonbonded/Sticky.cpp 2013/01/15 16:28:22 1833 @@ -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 @@ -45,64 +46,14 @@ #include #include "nonbonded/Sticky.hpp" #include "nonbonded/LJ.hpp" +#include "types/StickyAdapter.hpp" #include "utils/simError.h" using namespace std; namespace OpenMD { - - bool Sticky::initialized_ = false; - ForceField* Sticky::forceField_ = NULL; - map Sticky::StickyMap; - map, StickyInteractionData> Sticky::MixingMap; - Sticky* Sticky::_instance = NULL; - - Sticky* Sticky::Instance() { - if (!_instance) { - _instance = new Sticky(); - } - return _instance; - } - - StickyParam Sticky::getStickyParam(AtomType* atomType) { + Sticky::Sticky() : name_("Sticky"), initialized_(false), forceField_(NULL) {} - // Do sanity checking on the AtomType we were passed before - // building any data structures: - if (!atomType->isSticky() && !atomType->isStickyPower()) { - sprintf( painCave.errMsg, - "Sticky::getStickyParam was passed an atomType (%s) that does\n" - "\tnot appear to be a Sticky atom.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DirectionalAtomType* daType = dynamic_cast(atomType); - GenericData* data = daType->getPropertyByName("Sticky"); - if (data == NULL) { - sprintf( painCave.errMsg, "Sticky::getStickyParam could not find\n" - "\tSticky parameters for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - StickyParamGenericData* stickyData = dynamic_cast(data); - if (stickyData == NULL) { - sprintf( painCave.errMsg, - "Sticky::getStickyParam could not convert GenericData to\n" - "\tStickyParamGenericData for atom type %s\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - return stickyData->getData(); - } - void Sticky::initialize() { ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); @@ -115,8 +66,8 @@ namespace OpenMD { for (at = atomTypes->beginType(i); at != NULL; at = atomTypes->nextType(i)) { - if (at->isSticky() || at->isStickyPower()) - addType(at); + StickyAdapter sa = StickyAdapter(at); + if (sa.isSticky()) addType(at); } initialized_ = true; @@ -124,22 +75,20 @@ namespace OpenMD { void Sticky::addType(AtomType* atomType){ // add it to the map: - AtomTypeProperties atp = atomType->getATP(); pair::iterator,bool> ret; - ret = StickyMap.insert( pair(atp.ident, atomType) ); + ret = StickyMap.insert( pair(atomType->getIdent(), + atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "Sticky already had a previous entry with ident %d\n", - atp.ident); + atomType->getIdent() ); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); - } + } - RealType w0i, v0i, v0pi, rli, rui, rlpi, rupi; - - StickyParam sticky1 = getStickyParam(atomType); + StickyAdapter sticky1 = StickyAdapter(atomType); // Now, iterate over all known types and add to the mixing map: @@ -147,8 +96,8 @@ namespace OpenMD { for( it = StickyMap.begin(); it != StickyMap.end(); ++it) { AtomType* atype2 = (*it).second; - - StickyParam sticky2 = getStickyParam(atype2); + + StickyAdapter sticky2 = StickyAdapter(atype2); StickyInteractionData mixer; @@ -157,15 +106,15 @@ namespace OpenMD { // Lorentz- Berthelot mixing rules (which happen to do the right thing // when atomType and atype2 happen to be the same. - mixer.rl = 0.5 * ( sticky1.rl + sticky2.rl ); - mixer.ru = 0.5 * ( sticky1.ru + sticky2.ru ); - mixer.rlp = 0.5 * ( sticky1.rlp + sticky2.rlp ); - mixer.rup = 0.5 * ( sticky1.rup + sticky2.rup ); + mixer.rl = 0.5 * ( sticky1.getRl() + sticky2.getRl() ); + mixer.ru = 0.5 * ( sticky1.getRu() + sticky2.getRu() ); + mixer.rlp = 0.5 * ( sticky1.getRlp() + sticky2.getRlp() ); + mixer.rup = 0.5 * ( sticky1.getRup() + sticky2.getRup() ); mixer.rbig = max(mixer.ru, mixer.rup); - mixer.w0 = sqrt( sticky1.w0 * sticky2.w0 ); - mixer.v0 = sqrt( sticky1.v0 * sticky2.v0 ); - mixer.v0p = sqrt( sticky1.v0p * sticky2.v0p ); - mixer.isPower = atomType->isStickyPower() && atype2->isStickyPower(); + mixer.w0 = sqrt( sticky1.getW0() * sticky2.getW0() ); + mixer.v0 = sqrt( sticky1.getV0() * sticky2.getV0() ); + mixer.v0p = sqrt( sticky1.getV0p() * sticky2.getV0p() ); + mixer.isPower = sticky1.isStickyPower() && sticky2.isStickyPower(); CubicSpline* s = new CubicSpline(); s->addPoint(mixer.rl, 1.0); @@ -189,275 +138,217 @@ namespace OpenMD { } } - RealType Sticky::getStickyCut(int atid) { + /** + * This function does the sticky portion of the SSD potential + * [Chandra and Ichiye, Journal of Chemical Physics 111, 2701 + * (1999)]. The Lennard-Jones and dipolar interaction must be + * handled separately. We assume that the rotation matrices have + * already been calculated and placed in the A1 & A2 entries in the + * idat structure. + */ + + void Sticky::calcForce(InteractionData &idat) { + if (!initialized_) initialize(); - std::map :: const_iterator it; - it = StickyMap.find(atid); - if (it == StickyMap.end()) { - sprintf( painCave.errMsg, - "Sticky::getStickyCut could not find atid %d in StickyMap\n", - (atid)); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - AtomType* atype = it->second; - return MixingMap[make_pair(atype, atype)].rbig; - } - - - void Sticky::calcForce(AtomType* at1, AtomType* at2, Vector3d d, - RealType rij, RealType r2, RealType sw, - RealType &vpair, RealType &pot, - RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1, - Vector3d &t1, Vector3d &t2) { - - // This routine does only the sticky portion of the SSD potential - // [Chandra and Ichiye, J. Chem. Phys. 111, 2701 (1999)]. - // The Lennard-Jones and dipolar interaction must be handled separately. - // We assume that the rotation matrices have already been calculated - // and placed in the A array. - - if (!initialized_) initialize(); - - pair key = make_pair(at1, at2); - StickyInteractionData mixer = MixingMap[key]; + map, StickyInteractionData>::iterator it; + it = MixingMap.find(idat.atypes); + if (it != MixingMap.end()) { - RealType w0 = mixer.w0; - RealType v0 = mixer.v0; - RealType v0p = mixer.v0p; - RealType rl = mixer.rl; - RealType ru = mixer.ru; - RealType rlp = mixer.rlp; - RealType rup = mixer.rup; - RealType rbig = mixer.rbig; - bool isPower = mixer.isPower; - - if (rij <= rbig) { - - RealType r3 = r2 * rij; - RealType r5 = r3 * r2; - - RotMat3x3d A1trans = A1.transpose(); - RotMat3x3d A2trans = A2.transpose(); - - // rotate the inter-particle separation into the two different - // body-fixed coordinate systems: - - Vector3d ri = A1 * d; - - // negative sign because this is the vector from j to i: - - Vector3d rj = -A2 * d; - - RealType xi = ri.x(); - RealType yi = ri.y(); - RealType zi = ri.z(); - - RealType xj = rj.x(); - RealType yj = rj.y(); - RealType zj = rj.z(); - - RealType xi2 = xi * xi; - RealType yi2 = yi * yi; - RealType zi2 = zi * zi; - - RealType xj2 = xj * xj; - RealType yj2 = yj * yj; - RealType zj2 = zj * zj; - - // calculate the switching info. from the splines - - RealType s = 0.0; - RealType dsdr = 0.0; - RealType sp = 0.0; - RealType dspdr = 0.0; - - if (rij < ru) { - if (rij < rl) { - s = 1.0; - dsdr = 0.0; - } else { - // we are in the switching region - - pair res = mixer.s->getValueAndDerivativeAt(rij); - s = res.first; - dsdr = res.second; - } - } - - if (rij < rup) { - if (rij < rlp) { - sp = 1.0; - dspdr = 0.0; - } else { - // we are in the switching region - - pair res =mixer.sp->getValueAndDerivativeAt(rij); - sp = res.first; - dspdr = res.second; - } - } - - RealType wi = 2.0*(xi2-yi2)*zi / r3; - RealType wj = 2.0*(xj2-yj2)*zj / r3; - RealType w = wi+wj; - - - RealType zif = zi/rij - 0.6; - RealType zis = zi/rij + 0.8; - - RealType zjf = zj/rij - 0.6; - RealType zjs = zj/rij + 0.8; - - RealType wip = zif*zif*zis*zis - w0; - RealType wjp = zjf*zjf*zjs*zjs - w0; - RealType wp = wip + wjp; - - Vector3d dwi(4.0*xi*zi/r3 - 6.0*xi*zi*(xi2-yi2)/r5, - - 4.0*yi*zi/r3 - 6.0*yi*zi*(xi2-yi2)/r5, - 2.0*(xi2-yi2)/r3 - 6.0*zi2*(xi2-yi2)/r5); + StickyInteractionData mixer = (*it).second; - Vector3d dwj(4.0*xj*zj/r3 - 6.0*xj*zj*(xj2-yj2)/r5, - - 4.0*yj*zj/r3 - 6.0*yj*zj*(xj2-yj2)/r5, - 2.0*(xj2-yj2)/r3 - 6.0*zj2*(xj2-yj2)/r5); - - RealType uglyi = zif*zif*zis + zif*zis*zis; - RealType uglyj = zjf*zjf*zjs + zjf*zjs*zjs; - - Vector3d dwip(-2.0*xi*zi*uglyi/r3, - -2.0*yi*zi*uglyi/r3, - 2.0*(1.0/rij - zi2/r3)*uglyi); - - Vector3d dwjp(-2.0*xj*zj*uglyj/r3, - -2.0*yj*zj*uglyj/r3, - 2.0*(1.0/rij - zj2/r3)*uglyj); - - Vector3d dwidu(4.0*(yi*zi2 + 0.5*yi*(xi2-yi2))/r3, - 4.0*(xi*zi2 - 0.5*xi*(xi2-yi2))/r3, - - 8.0*xi*yi*zi/r3); - - Vector3d dwjdu(4.0*(yj*zj2 + 0.5*yj*(xj2-yj2))/r3, - 4.0*(xj*zj2 - 0.5*xj*(xj2-yj2))/r3, - - 8.0*xj*yj*zj/r3); + RealType w0 = mixer.w0; + RealType v0 = mixer.v0; + RealType v0p = mixer.v0p; + RealType rl = mixer.rl; + RealType ru = mixer.ru; + RealType rlp = mixer.rlp; + RealType rup = mixer.rup; + RealType rbig = mixer.rbig; + bool isPower = mixer.isPower; - Vector3d dwipdu(2.0*yi*uglyi/rij, - -2.0*xi*uglyi/rij, - 0.0); - - Vector3d dwjpdu(2.0*yj*uglyj/rij, - -2.0*xj*uglyj/rij, - 0.0); + if ( *(idat.rij) <= rbig) { + + RealType r3 = *(idat.r2) * *(idat.rij); + RealType r5 = r3 * *(idat.r2); + + RotMat3x3d A1trans = idat.A1->transpose(); + RotMat3x3d A2trans = idat.A2->transpose(); + + // rotate the inter-particle separation into the two different + // body-fixed coordinate systems: + + Vector3d ri = *(idat.A1) * *(idat.d); + + // negative sign because this is the vector from j to i: + + Vector3d rj = - *(idat.A2) * *(idat.d); + + RealType xi = ri.x(); + RealType yi = ri.y(); + RealType zi = ri.z(); + + RealType xj = rj.x(); + RealType yj = rj.y(); + RealType zj = rj.z(); + + RealType xi2 = xi * xi; + RealType yi2 = yi * yi; + RealType zi2 = zi * zi; + + RealType xj2 = xj * xj; + RealType yj2 = yj * yj; + RealType zj2 = zj * zj; + + // calculate the switching info. from the splines + + RealType s = 0.0; + RealType dsdr = 0.0; + RealType sp = 0.0; + RealType dspdr = 0.0; + + if ( *(idat.rij) < ru) { + if ( *(idat.rij) < rl) { + s = 1.0; + dsdr = 0.0; + } else { + // we are in the switching region + + pair res = mixer.s->getValueAndDerivativeAt(*(idat.rij)); + s = res.first; + dsdr = res.second; + } + } + + if (*(idat.rij) < rup) { + if ( *(idat.rij) < rlp) { + sp = 1.0; + dspdr = 0.0; + } else { + // we are in the switching region + + pair res =mixer.sp->getValueAndDerivativeAt( *(idat.rij)); + sp = res.first; + dspdr = res.second; + } + } + + RealType wi = 2.0*(xi2-yi2)*zi / r3; + RealType wj = 2.0*(xj2-yj2)*zj / r3; + RealType w = wi+wj; + + + RealType zif = zi/ *(idat.rij) - 0.6; + RealType zis = zi/ *(idat.rij) + 0.8; + + RealType zjf = zj/ *(idat.rij) - 0.6; + RealType zjs = zj/ *(idat.rij) + 0.8; + + RealType wip = zif*zif*zis*zis - w0; + RealType wjp = zjf*zjf*zjs*zjs - w0; + RealType wp = wip + wjp; + + Vector3d dwi(4.0*xi*zi/r3 - 6.0*xi*zi*(xi2-yi2)/r5, + - 4.0*yi*zi/r3 - 6.0*yi*zi*(xi2-yi2)/r5, + 2.0*(xi2-yi2)/r3 - 6.0*zi2*(xi2-yi2)/r5); + + Vector3d dwj(4.0*xj*zj/r3 - 6.0*xj*zj*(xj2-yj2)/r5, + - 4.0*yj*zj/r3 - 6.0*yj*zj*(xj2-yj2)/r5, + 2.0*(xj2-yj2)/r3 - 6.0*zj2*(xj2-yj2)/r5); + + RealType uglyi = zif*zif*zis + zif*zis*zis; + RealType uglyj = zjf*zjf*zjs + zjf*zjs*zjs; - if (isPower) { - RealType frac1 = 0.25; - RealType frac2 = 0.75; - RealType wi2 = wi*wi; - RealType wj2 = wj*wj; - // sticky power has no w' function: - w = frac1 * wi * wi2 + frac2*wi + frac1*wj*wj2 + frac2*wj + v0p; - wp = 0.0; - dwi = frac1*3.0*wi2*dwi + frac2*dwi; - dwj = frac1*3.0*wj2*dwi + frac2*dwi; - dwip = V3Zero; - dwjp = V3Zero; - dwidu = frac1*3.0*wi2*dwidu + frac2*dwidu; - dwidu = frac1*3.0*wj2*dwjdu + frac2*dwjdu; - dwipdu = V3Zero; - dwjpdu = V3Zero; - sp = 0.0; - dspdr = 0.0; - } + Vector3d dwip(-2.0*xi*zi*uglyi/r3, + -2.0*yi*zi*uglyi/r3, + 2.0*(1.0/ *(idat.rij) - zi2/r3)*uglyi); + + Vector3d dwjp(-2.0*xj*zj*uglyj/r3, + -2.0*yj*zj*uglyj/r3, + 2.0*(1.0/ *(idat.rij) - zj2/r3)*uglyj); + + Vector3d dwidu(4.0*(yi*zi2 + 0.5*yi*(xi2-yi2))/r3, + 4.0*(xi*zi2 - 0.5*xi*(xi2-yi2))/r3, + - 8.0*xi*yi*zi/r3); + + Vector3d dwjdu(4.0*(yj*zj2 + 0.5*yj*(xj2-yj2))/r3, + 4.0*(xj*zj2 - 0.5*xj*(xj2-yj2))/r3, + - 8.0*xj*yj*zj/r3); + + Vector3d dwipdu(2.0*yi*uglyi/ *(idat.rij) , + -2.0*xi*uglyi/ *(idat.rij) , + 0.0); + + Vector3d dwjpdu(2.0*yj*uglyj/ *(idat.rij) , + -2.0*xj*uglyj/ *(idat.rij) , + 0.0); + + if (isPower) { + RealType frac1 = 0.25; + RealType frac2 = 0.75; + RealType wi2 = wi*wi; + RealType wj2 = wj*wj; + // sticky power has no w' function: + w = frac1 * wi * wi2 + frac2*wi + frac1*wj*wj2 + frac2*wj + v0p; + wp = 0.0; + dwi = frac1*RealType(3.0)*wi2*dwi + frac2*dwi; + dwj = frac1*RealType(3.0)*wj2*dwi + frac2*dwi; + dwip = V3Zero; + dwjp = V3Zero; + dwidu = frac1*RealType(3.0)*wi2*dwidu + frac2*dwidu; + dwidu = frac1*RealType(3.0)*wj2*dwjdu + frac2*dwjdu; + dwipdu = V3Zero; + dwjpdu = V3Zero; + sp = 0.0; + dspdr = 0.0; + } + - vpair += 0.5*(v0*s*w + v0p*sp*wp); - pot += 0.5*(v0*s*w + v0p*sp*wp)*sw; - // do the torques first since they are easy: - // remember that these are still in the body-fixed axes - - Vector3d ti = 0.5*sw*(v0*s*dwidu + v0p*sp*dwipdu); - Vector3d tj = 0.5*sw*(v0*s*dwjdu + v0p*sp*dwjpdu); - - // go back to lab frame using transpose of rotation matrix: - - t1 += A1trans * ti; - t2 += A2trans * tj; - - // Now, on to the forces: - - // first rotate the i terms back into the lab frame: - - Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) * sw; - Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) * sw; - - Vector3d fii = A1trans * radcomi; - Vector3d fjj = A2trans * radcomj; - - // now assemble these with the radial-only terms: - - f1 += 0.5 * ((v0*dsdr*w + v0p*dspdr*wp) * d / rij + fii - fjj); - + *(idat.vpair) += RealType(0.5)*(v0*s*w + v0p*sp*wp); + (*(idat.pot))[HYDROGENBONDING_FAMILY] += RealType(0.5)*(v0*s*w + v0p*sp*wp)* *(idat.sw) ; + + // do the torques first since they are easy: + // remember that these are still in the body-fixed axes + + Vector3d ti = RealType(0.5)* *(idat.sw) *(v0*s*dwidu + v0p*sp*dwipdu); + Vector3d tj = RealType(0.5)* *(idat.sw) *(v0*s*dwjdu + v0p*sp*dwjpdu); + + // go back to lab frame using transpose of rotation matrix: + + *(idat.t1) += A1trans * ti; + *(idat.t2) += A2trans * tj; + + // Now, on to the forces: + + // first rotate the i terms back into the lab frame: + + Vector3d radcomi = (v0 * s * dwi + v0p * sp * dwip) * *(idat.sw); + Vector3d radcomj = (v0 * s * dwj + v0p * sp * dwjp) * *(idat.sw); + + Vector3d fii = A1trans * radcomi; + Vector3d fjj = A2trans * radcomj; + + // now assemble these with the radial-only terms: + + *(idat.f1) += RealType(0.5) * ((v0*dsdr*w + v0p*dspdr*wp) * *(idat.d) / + *(idat.rij) + fii - fjj); + + } } - - return; + return; } - void Sticky::do_sticky_pair(int *atid1, int *atid2, RealType *d, - RealType *r, RealType *r2, RealType *sw, - RealType *vpair, RealType *pot, RealType *A1, - RealType *A2, RealType *f1, - RealType *t1, RealType *t2) { - - if (!initialized_) initialize(); - - AtomType* atype1 = StickyMap[*atid1]; - AtomType* atype2 = StickyMap[*atid2]; - - Vector3d disp(d); - Vector3d frc(f1); - Vector3d trq1(t1); - Vector3d trq2(t2); - RotMat3x3d Ai(A1); - RotMat3x3d Aj(A2); - - calcForce(atype1, atype2, disp, *r, *r2, *sw, *vpair, *pot, - Ai, Aj, frc, trq1, trq2); - - f1[0] = frc.x(); - f1[1] = frc.y(); - f1[2] = frc.z(); - - t1[0] = trq1.x(); - t1[1] = trq1.y(); - t1[2] = trq1.z(); - - t2[0] = trq2.x(); - t2[1] = trq2.y(); - t2[2] = trq2.z(); - - return; + RealType Sticky::getSuggestedCutoffRadius(pair atypes) { + if (!initialized_) initialize(); + map, StickyInteractionData>::iterator it; + it = MixingMap.find(atypes); + if (it == MixingMap.end()) + return 0.0; + else { + StickyInteractionData mixer = (*it).second; + return mixer.rbig; + } } } - -extern "C" { - -#define fortranGetStickyCut FC_FUNC(getstickycut, GETSTICKYCUT) -#define fortranDoStickyPair FC_FUNC(do_sticky_pair, DO_STICKY_PAIR) - - RealType fortranGetStickyCut(int* atid) { - return OpenMD::Sticky::Instance()->getStickyCut(*atid); - } - - void fortranDoStickyPair(int *atid1, int *atid2, RealType *d, RealType *r, - RealType *r2, RealType *sw, RealType *vpair, RealType *pot, - RealType *A1, RealType *A2, RealType *f1, - RealType *t1, RealType *t2){ - - return OpenMD::Sticky::Instance()->do_sticky_pair(atid1, atid2, d, r, r2, - sw, vpair, pot, A1, A2, - f1, t1, t2); - } -}