--- branches/development/src/nonbonded/MAW.cpp 2011/01/05 14:49:05 1536 +++ branches/development/src/nonbonded/MAW.cpp 2011/06/16 22:00:08 1583 @@ -50,8 +50,7 @@ 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() { @@ -128,9 +127,8 @@ namespace OpenMD { 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 +143,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 +171,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 +196,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 +215,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[0] += pot_temp; - idat.pot[0] += 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 +235,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 +264,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 {