--- branches/development/src/nonbonded/MAW.cpp 2010/12/29 19:59:21 1532 +++ branches/development/src/nonbonded/MAW.cpp 2011/11/22 20:38:56 1665 @@ -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 @@ -50,53 +51,43 @@ namespace OpenMD { namespace OpenMD { - MAW::MAW() : name_("MAW"), initialized_(false), forceField_(NULL), - shiftedPot_(false), shiftedFrc_(false) {} + MAW::MAW() : name_("MAW"), initialized_(false), forceField_(NULL) {} void MAW::initialize() { ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes(); ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j; NonBondedInteractionType* nbt; + ForceField::NonBondedInteractionTypeContainer::KeyType keys; for (nbt = nbiTypes->beginType(j); nbt != NULL; nbt = nbiTypes->nextType(j)) { if (nbt->isMAW()) { - pair atypes = nbt->getAtomTypes(); - - GenericData* data = nbt->getPropertyByName("MAW"); - if (data == NULL) { - sprintf( painCave.errMsg, "MAW::initialize could not find\n" - "\tMAW parameters for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - MAWData* mawData = dynamic_cast(data); - if (mawData == NULL) { + keys = nbiTypes->getKeys(j); + AtomType* at1 = forceField_->getAtomType(keys[0]); + AtomType* at2 = forceField_->getAtomType(keys[1]); + + MAWInteractionType* mit = dynamic_cast(nbt); + + if (mit == NULL) { sprintf( painCave.errMsg, - "MAW::initialize could not convert GenericData to\n" - "\tMAWData for %s - %s interaction.\n", - atypes.first->getName().c_str(), - atypes.second->getName().c_str()); + "MAW::initialize could not convert NonBondedInteractionType\n" + "\tto MAWInteractionType for %s - %s interaction.\n", + at1->getName().c_str(), + at2->getName().c_str()); painCave.severity = OPENMD_ERROR; painCave.isFatal = 1; simError(); } - MAWParam mawParam = mawData->getData(); - - RealType De = mawParam.De; - RealType beta = mawParam.beta; - RealType Re = mawParam.Re; - RealType ca1 = mawParam.ca1; - RealType cb1 = mawParam.cb1; + RealType De = mit->getD(); + RealType beta = mit->getBeta(); + RealType Re = mit->getR(); + RealType ca1 = mit->getCA1(); + RealType cb1 = mit->getCB1(); - addExplicitInteraction(atypes.first, atypes.second, + addExplicitInteraction(at1, at2, De, beta, Re, ca1, cb1); } } @@ -124,13 +115,12 @@ namespace OpenMD { } } - void MAW::calcForce(InteractionData idat) { + void MAW::calcForce(InteractionData &idat) { if (!initialized_) initialize(); - pair key = make_pair(idat.atype1, idat.atype2); map, MAWInteractionData>::iterator it; - it = MixingMap.find(key); + it = MixingMap.find( idat.atypes ); if (it != MixingMap.end()) { MAWInteractionData mixer = (*it).second; @@ -145,24 +135,24 @@ namespace OpenMD { RealType ca1 = mixer.ca1; RealType cb1 = mixer.cb1; - bool j_is_Metal = idat.atype2->isMetal(); + bool j_is_Metal = idat.atypes.second->isMetal(); Vector3d r; RotMat3x3d Atrans; if (j_is_Metal) { // rotate the inter-particle separation into the two different // body-fixed coordinate systems: - r = idat.A1 * idat.d; - Atrans = idat.A1.transpose(); + r = *(idat.A1) * *(idat.d); + Atrans = idat.A1->transpose(); } else { // negative sign because this is the vector from j to i: - r = -idat.A2 * idat.d; - Atrans = idat.A2.transpose(); + r = -*(idat.A2) * *(idat.d); + Atrans = idat.A2->transpose(); } // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2) - RealType expt = -beta*(idat.rij - R_e); + RealType expt = -beta*( *(idat.rij) - R_e); RealType expfnc = exp(expt); RealType expfnc2 = expfnc*expfnc; @@ -173,19 +163,19 @@ namespace OpenMD { myPot = D_e * (expfnc2 - 2.0 * expfnc); myDeriv = 2.0 * D_e * beta * (expfnc - expfnc2); - if (MAW::shiftedPot_ || MAW::shiftedFrc_) { - exptC = -beta*(idat.rcut - R_e); + if (idat.shiftedPot || idat.shiftedForce) { + exptC = -beta*( *(idat.rcut) - R_e); expfncC = exp(exptC); expfnc2C = expfncC*expfncC; } - if (MAW::shiftedPot_) { + if (idat.shiftedPot) { myPotC = D_e * (expfnc2C - 2.0 * expfncC); myDerivC = 0.0; - } else if (MAW::shiftedFrc_) { + } else if (idat.shiftedForce) { myPotC = D_e * (expfnc2C - 2.0 * expfncC); myDerivC = 2.0 * D_e * beta * (expfnc2C - expfnc2C); - myPotC += myDerivC * (idat.rij - idat.rcut); + myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) ); } else { myPotC = 0.0; myDerivC = 0.0; @@ -198,8 +188,8 @@ namespace OpenMD { RealType y2 = y * y; RealType z2 = z * z; - RealType r3 = idat.r2 * idat.rij; - RealType r4 = idat.r2 * idat.r2; + RealType r3 = *(idat.r2) * *(idat.rij) ; + RealType r4 = *(idat.r2) * *(idat.r2); // angular modulation of morse part of potential to approximate // the squares of the two HOMO lone pair orbitals in water: @@ -217,18 +207,19 @@ namespace OpenMD { // Vmorse(r)*[a*x2/r2 + b*z/r + (1-a-b)] RealType Vmorse = (myPot - myPotC); - RealType Vang = ca1 * x2 / idat.r2 + cb1 * z / idat.rij + (0.8 - ca1 / 3.0); - - RealType pot_temp = idat.vdwMult * Vmorse * Vang; - idat.vpair += pot_temp; - idat.pot += idat.sw * pot_temp; + RealType Vang = ca1 * x2 / *(idat.r2) + + cb1 * z / *(idat.rij) + (0.8 - ca1 / 3.0); + + RealType pot_temp = *(idat.vdwMult) * Vmorse * Vang; + *(idat.vpair) += pot_temp; + (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp; - Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) / idat.rij; - - Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / idat.r2 - cb1 * x * z / r3, - -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3, - -2.0 * ca1 * x2 * z / r4 + cb1 / idat.rij - cb1 * z2 / r3); + Vector3d dVmorsedr = (myDeriv - myDerivC) * Vector3d(x, y, z) / *(idat.rij) ; + Vector3d dVangdr = Vector3d(-2.0 * ca1 * x2 * x / r4 + 2.0 * ca1 * x / *(idat.r2) - cb1 * x * z / r3, + -2.0 * ca1 * x2 * y / r4 - cb1 * z * y / r3, + -2.0 * ca1 * x2 * z / r4 + cb1 / *(idat.rij) - cb1 * z2 / r3); + // chain rule to put these back on x, y, z Vector3d dvdr = Vang * dVmorsedr + Vmorse * dVangdr; @@ -236,26 +227,26 @@ namespace OpenMD { // Torques for Vang using method of Price: // S. L. Price, A. J. Stone, and M. Alderton, Mol. Phys. 52, 987 (1984). - Vector3d dVangdu = Vector3d(cb1 * y / idat.rij, - 2.0 * ca1 * x * z / idat.r2 - cb1 * x / idat.rij, - -2.0 * ca1 * y * x / idat.r2); + Vector3d dVangdu = Vector3d(cb1 * y / *(idat.rij) , + 2.0 * ca1 * x * z / *(idat.r2) - cb1 * x / *(idat.rij), + -2.0 * ca1 * y * x / *(idat.r2)); // do the torques first since they are easy: // remember that these are still in the body fixed axes - Vector3d trq = idat.vdwMult * Vmorse * dVangdu * idat.sw; + Vector3d trq = *(idat.vdwMult) * Vmorse * dVangdu * *(idat.sw); // go back to lab frame using transpose of rotation matrix: if (j_is_Metal) { - idat.t1 += Atrans * trq; + *(idat.t1) += Atrans * trq; } else { - idat.t2 += Atrans * trq; + *(idat.t2) += Atrans * trq; } // Now, on to the forces (still in body frame of water) - Vector3d ftmp = idat.vdwMult * idat.sw * dvdr; + Vector3d ftmp = *(idat.vdwMult) * *(idat.sw) * dvdr; // rotate the terms back into the lab frame: Vector3d flab; @@ -265,17 +256,16 @@ namespace OpenMD { flab = - Atrans * ftmp; } - idat.f1 += flab; + *(idat.f1) += flab; } return; } - RealType MAW::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) { + RealType MAW::getSuggestedCutoffRadius(pair atypes) { if (!initialized_) initialize(); - pair key = make_pair(at1, at2); map, MAWInteractionData>::iterator it; - it = MixingMap.find(key); + it = MixingMap.find(atypes); if (it == MixingMap.end()) return 0.0; else {