--- branches/development/src/nonbonded/Electrostatic.cpp 2011/11/22 20:38:56 1665 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/05/18 21:44:02 1710 @@ -47,7 +47,8 @@ #include "nonbonded/Electrostatic.hpp" #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" -#include "types/DirectionalAtomType.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/MultipoleAdapter.hpp" #include "io/Globals.hpp" namespace OpenMD { @@ -280,132 +281,42 @@ namespace OpenMD { electrostaticAtomData.is_SplitDipole = false; electrostaticAtomData.is_Quadrupole = false; - if (atomType->isCharge()) { - GenericData* data = atomType->getPropertyByName("Charge"); + FixedChargeAdapter fca = FixedChargeAdapter(atomType); - if (data == NULL) { - sprintf( painCave.errMsg, "Electrostatic::addType could not find " - "Charge\n" - "\tparameters for atomType %s.\n", - atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DoubleGenericData* doubleData = dynamic_cast(data); - if (doubleData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "Charge for\n" - "\tatom type %s\n", atomType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + if (fca.isFixedCharge()) { electrostaticAtomData.is_Charge = true; - electrostaticAtomData.charge = doubleData->getData(); + electrostaticAtomData.charge = fca.getCharge(); } - if (atomType->isDirectional()) { - DirectionalAtomType* daType = dynamic_cast(atomType); - - if (daType->isDipole()) { - GenericData* data = daType->getPropertyByName("Dipole"); - - if (data == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not find Dipole\n" - "\tparameters for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DoubleGenericData* doubleData = dynamic_cast(data); - if (doubleData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "Dipole Moment\n" - "\tfor atom type %s\n", daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + MultipoleAdapter ma = MultipoleAdapter(atomType); + if (ma.isMultipole()) { + if (ma.isDipole()) { electrostaticAtomData.is_Dipole = true; - electrostaticAtomData.dipole_moment = doubleData->getData(); + electrostaticAtomData.dipole_moment = ma.getDipoleMoment(); } - - if (daType->isSplitDipole()) { - GenericData* data = daType->getPropertyByName("SplitDipoleDistance"); - - if (data == NULL) { - sprintf(painCave.errMsg, - "Electrostatic::addType could not find SplitDipoleDistance\n" - "\tparameter for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - - DoubleGenericData* doubleData = dynamic_cast(data); - if (doubleData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "SplitDipoleDistance for\n" - "\tatom type %s\n", daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } + if (ma.isSplitDipole()) { electrostaticAtomData.is_SplitDipole = true; - electrostaticAtomData.split_dipole_distance = doubleData->getData(); + electrostaticAtomData.split_dipole_distance = ma.getSplitDipoleDistance(); } - - if (daType->isQuadrupole()) { - GenericData* data = daType->getPropertyByName("QuadrupoleMoments"); - - if (data == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not find QuadrupoleMoments\n" - "\tparameter for atomType %s.\n", - daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - + if (ma.isQuadrupole()) { // Quadrupoles in OpenMD are set as the diagonal elements // of the diagonalized traceless quadrupole moment tensor. // The column vectors of the unitary matrix that diagonalizes // the quadrupole moment tensor become the eFrame (or the // electrostatic version of the body-fixed frame. - - Vector3dGenericData* v3dData = dynamic_cast(data); - if (v3dData == NULL) { - sprintf( painCave.errMsg, - "Electrostatic::addType could not convert GenericData to " - "Quadrupole Moments for\n" - "\tatom type %s\n", daType->getName().c_str()); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } electrostaticAtomData.is_Quadrupole = true; - electrostaticAtomData.quadrupole_moments = v3dData->getData(); + electrostaticAtomData.quadrupole_moments = ma.getQuadrupoleMoments(); } } - AtomTypeProperties atp = atomType->getATP(); pair::iterator,bool> ret; - ret = ElectrostaticList.insert( pair(atp.ident, atomType) ); + ret = ElectrostaticList.insert( pair(atomType->getIdent(), + atomType) ); if (ret.second == false) { sprintf( painCave.errMsg, "Electrostatic already had a previous entry with ident %d\n", - atp.ident); + atomType->getIdent() ); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); @@ -461,6 +372,7 @@ namespace OpenMD { RealType c1, c2, c3, c4; RealType erfcVal(1.0), derfcVal(0.0); RealType BigR; + RealType two(2.0), three(3.0); Vector3d Q_i, Q_j; Vector3d ux_i, uy_i, uz_i; @@ -621,7 +533,7 @@ namespace OpenMD { if (idat.excluded) { indirect_vpair += preVal * rfVal; indirect_Pot += *(idat.sw) * preVal * rfVal; - indirect_dVdr += *(idat.sw) * preVal * 2.0 * rfVal * riji * rhat; + indirect_dVdr += *(idat.sw) * preVal * two * rfVal * riji * rhat; } } else { @@ -649,7 +561,7 @@ namespace OpenMD { vpair += vterm; epot += *(idat.sw) * vterm; - dVdr += -preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); + dVdr += -preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j); duduz_j += -preSw * rhat * (ri2 - preRF2_ * *(idat.rij) ); // Even if we excluded this pair from direct interactions, @@ -738,7 +650,7 @@ namespace OpenMD { c2ri = c2 * riji; c3ri = c3 * riji; c4rij = c4 * *(idat.rij) ; - rhatdot2 = 2.0 * rhat * c3; + rhatdot2 = two * rhat * c3; rhatc4 = rhat * c4rij; // calculate the potential @@ -751,9 +663,9 @@ namespace OpenMD { // calculate derivatives for the forces and torques - dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (2.0*cx_j*ux_j + rhat)*c3ri) + - qyy_j* (cy2*rhatc4 - (2.0*cy_j*uy_j + rhat)*c3ri) + - qzz_j* (cz2*rhatc4 - (2.0*cz_j*uz_j + rhat)*c3ri)); + dVdr += -preSw * ( qxx_j* (cx2*rhatc4 - (two*cx_j*ux_j + rhat)*c3ri) + + qyy_j* (cy2*rhatc4 - (two*cy_j*uy_j + rhat)*c3ri) + + qzz_j* (cz2*rhatc4 - (two*cz_j*uz_j + rhat)*c3ri)); dudux_j += preSw * qxx_j * cx_j * rhatdot2; duduy_j += preSw * qyy_j * cy_j * rhatdot2; @@ -777,7 +689,7 @@ namespace OpenMD { vpair += vterm; epot += *(idat.sw) * vterm; - dVdr += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_ * uz_i); + dVdr += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_ * uz_i); duduz_i += preSw * rhat * (ri2 - preRF2_ * *(idat.rij) ); @@ -855,10 +767,10 @@ namespace OpenMD { a1 = 5.0 * ct_i * ct_j - ct_ij; - dVdr += preSw * 3.0 * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i); + dVdr += preSw * three * ri4 * (a1 * rhat - ct_i * uz_j - ct_j * uz_i); - duduz_i += preSw * (ri3 * (uz_j - 3.0 * ct_j * rhat) - preRF2_*uz_j); - duduz_j += preSw * (ri3 * (uz_i - 3.0 * ct_i * rhat) - preRF2_*uz_i); + duduz_i += preSw * (ri3 * (uz_j - three * ct_j * rhat) - preRF2_*uz_j); + duduz_j += preSw * (ri3 * (uz_i - three * ct_i * rhat) - preRF2_*uz_i); if (idat.excluded) { indirect_vpair += - pref * preRF2_ * ct_ij; @@ -963,7 +875,7 @@ namespace OpenMD { c2ri = c2 * riji; c3ri = c3 * riji; c4rij = c4 * *(idat.rij) ; - rhatdot2 = 2.0 * rhat * c3; + rhatdot2 = two * rhat * c3; rhatc4 = rhat * c4rij; // calculate the potential @@ -977,9 +889,9 @@ namespace OpenMD { // calculate the derivatives for the forces and torques - dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (2.0*cx_i*ux_i + rhat)*c3ri) + - qyy_i* (cy2*rhatc4 - (2.0*cy_i*uy_i + rhat)*c3ri) + - qzz_i* (cz2*rhatc4 - (2.0*cz_i*uz_i + rhat)*c3ri)); + dVdr += -preSw * (qxx_i* (cx2*rhatc4 - (two*cx_i*ux_i + rhat)*c3ri) + + qyy_i* (cy2*rhatc4 - (two*cy_i*uy_i + rhat)*c3ri) + + qzz_i* (cz2*rhatc4 - (two*cz_i*uz_i + rhat)*c3ri)); dudux_i += preSw * qxx_i * cx_i * rhatdot2; duduy_i += preSw * qyy_i * cy_i * rhatdot2;