--- branches/development/src/nonbonded/Electrostatic.cpp 2012/01/06 19:03:05 1668 +++ branches/development/src/nonbonded/Electrostatic.cpp 2012/05/24 01:48:29 1720 @@ -47,9 +47,14 @@ #include "nonbonded/Electrostatic.hpp" #include "utils/simError.h" #include "types/NonBondedInteractionType.hpp" -#include "types/DirectionalAtomType.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/FluctuatingChargeAdapter.hpp" +#include "types/MultipoleAdapter.hpp" #include "io/Globals.hpp" +#include "nonbonded/SlaterIntegrals.hpp" +#include "utils/PhysicalConstants.hpp" + namespace OpenMD { Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), @@ -280,138 +285,102 @@ 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.fixedCharge = 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(); + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); + if (fqa.isFluctuatingCharge()) { + electrostaticAtomData.is_Fluctuating = true; + electrostaticAtomData.electronegativity = fqa.getElectronegativity(); + electrostaticAtomData.hardness = fqa.getHardness(); + electrostaticAtomData.slaterN = fqa.getSlaterN(); + electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); + } + 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(); } - ElectrostaticMap[atomType] = electrostaticAtomData; + ElectrostaticMap[atomType] = electrostaticAtomData; + + // Now, iterate over all known types and add to the mixing map: + + map::iterator it; + for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) { + AtomType* atype2 = (*it).first; + ElectrostaticAtomData eaData2 = (*it).second; + if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) { + + RealType a = electrostaticAtomData.slaterZeta; + RealType b = eaData2.slaterZeta; + int m = electrostaticAtomData.slaterN; + int n = eaData2.slaterN; + + // Create the spline of the coulombic integral for s-type + // Slater orbitals. Add a 2 angstrom safety window to deal + // with cutoffGroups that have charged atoms longer than the + // cutoffRadius away from each other. + + RealType rval; + RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1); + vector rvals; + vector J1vals; + vector J2vals; + for (int i = 0; i < np_; i++) { + rval = RealType(i) * dr; + rvals.push_back(rval); + J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); + J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); + } + + CubicSpline* J1 = new CubicSpline(); + J1->addPoints(rvals, J1vals); + CubicSpline* J2 = new CubicSpline(); + J2->addPoints(rvals, J2vals); + + pair key1, key2; + key1 = make_pair(atomType, atype2); + key2 = make_pair(atype2, atomType); + + Jij[key1] = J1; + Jij[key2] = J2; + } + } + return; } @@ -502,7 +471,7 @@ namespace OpenMD { bool j_is_Quadrupole = data2.is_Quadrupole; if (i_is_Charge) { - q_i = data1.charge; + q_i = data1.fixedCharge; if (idat.excluded) { *(idat.skippedCharge2) += q_i; } @@ -540,7 +509,7 @@ namespace OpenMD { } if (j_is_Charge) { - q_j = data2.charge; + q_j = data2.fixedCharge; if (idat.excluded) { *(idat.skippedCharge1) += q_j; } @@ -1053,7 +1022,7 @@ namespace OpenMD { } } else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) { if (i_is_Charge) { - chg1 = data.charge; + chg1 = data.fixedCharge; if (screeningMethod_ == DAMPED) { self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_; } else {