--- branches/development/src/nonbonded/Sticky.cpp 2011/11/22 20:38:56 1665 +++ branches/development/src/nonbonded/Sticky.cpp 2013/01/15 16:28:22 1833 @@ -46,52 +46,14 @@ #include #include "nonbonded/Sticky.hpp" #include "nonbonded/LJ.hpp" +#include "types/StickyAdapter.hpp" #include "utils/simError.h" using namespace std; namespace OpenMD { Sticky::Sticky() : name_("Sticky"), initialized_(false), forceField_(NULL) {} - - StickyParam Sticky::getStickyParam(AtomType* atomType) { - // 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(); @@ -104,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; @@ -113,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: @@ -136,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; @@ -146,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); @@ -331,26 +291,28 @@ namespace OpenMD { // 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; + dwi = frac1*RealType(3.0)*wi2*dwi + frac2*dwi; + dwj = frac1*RealType(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; + 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; } - *(idat.vpair) += 0.5*(v0*s*w + v0p*sp*wp); - (*(idat.pot))[HYDROGENBONDING_FAMILY] += 0.5*(v0*s*w + v0p*sp*wp)* *(idat.sw) ; + + + *(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 = 0.5* *(idat.sw) *(v0*s*dwidu + v0p*sp*dwipdu); - Vector3d tj = 0.5* *(idat.sw) *(v0*s*dwjdu + v0p*sp*dwjpdu); + 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: @@ -369,8 +331,8 @@ namespace OpenMD { // now assemble these with the radial-only terms: - *(idat.f1) += 0.5 * ((v0*dsdr*w + v0p*dspdr*wp) * *(idat.d) / - *(idat.rij) + fii - fjj); + *(idat.f1) += RealType(0.5) * ((v0*dsdr*w + v0p*dspdr*wp) * *(idat.d) / + *(idat.rij) + fii - fjj); } }