--- trunk/src/nonbonded/Electrostatic.cpp 2013/06/16 15:15:42 1879 +++ trunk/src/nonbonded/Electrostatic.cpp 2013/07/01 21:09:37 1895 @@ -213,16 +213,28 @@ namespace OpenMD { haveDampingAlpha_ = true; } - // find all of the Electrostatic atom Types: - ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); - ForceField::AtomTypeContainer::MapTypeIterator i; - AtomType* at; + Etypes.clear(); + Etids.clear(); + FQtypes.clear(); + FQtids.clear(); + ElectrostaticMap.clear(); + Jij.clear(); + nElectro_ = 0; + nFlucq_ = 0; + + Etids.resize( forceField_->getNAtomType(), -1); + FQtids.resize( forceField_->getNAtomType(), -1); + + set::iterator at; + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isElectrostatic()) nElectro_++; + if ((*at)->isFluctuatingCharge()) nFlucq_++; + } - for (at = atomTypes->beginType(i); at != NULL; - at = atomTypes->nextType(i)) { - - if (at->isElectrostatic()) - addType(at); + Jij.resize(nFlucq_); + + for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { + if ((*at)->isElectrostatic()) addType(*at); } if (summationMethod_ == esm_REACTION_FIELD) { @@ -398,7 +410,7 @@ namespace OpenMD { v11 = g - gc - rmRc*hc; v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric; v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric); - v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric; + v31 = (h-g*ri)*ri - (hc-gc*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric; v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric) - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric); v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2 @@ -455,7 +467,7 @@ namespace OpenMD { v11 = g - gc; v21 = g*ri - gc*ric; v22 = h - g*ri - (hc - gc*ric); - v31 = (h-g*ri)*ri - (hc-g*ric)*ric; + v31 = (h-g*ri)*ri - (hc-gc*ric)*ric; v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric); v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2; v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric; @@ -515,7 +527,7 @@ namespace OpenMD { v11 = g - gc; v21 = g*ri - gc*ric; v22 = h - g*ri - (hc - gc*ric); - v31 = (h-g*ri)*ri - (hc-g*ric)*ric; + v31 = (h-g*ri)*ri - (hc-gc*ric)*ric; v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric); v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2; v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric; @@ -625,7 +637,7 @@ namespace OpenMD { } void Electrostatic::addType(AtomType* atomType){ - + ElectrostaticAtomData electrostaticAtomData; electrostaticAtomData.is_Charge = false; electrostaticAtomData.is_Dipole = false; @@ -661,38 +673,54 @@ namespace OpenMD { electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); } - pair::iterator,bool> ret; - ret = ElectrostaticList.insert( pair(atomType->getIdent(), - atomType) ); + int atid = atomType->getIdent(); + int etid = Etypes.size(); + int fqtid = FQtypes.size(); + + pair::iterator,bool> ret; + ret = Etypes.insert( atid ); if (ret.second == false) { sprintf( painCave.errMsg, "Electrostatic already had a previous entry with ident %d\n", - atomType->getIdent() ); + atid); painCave.severity = OPENMD_INFO; painCave.isFatal = 0; simError(); } - ElectrostaticMap[atomType] = electrostaticAtomData; + Etids[ atid ] = etid; + ElectrostaticMap.push_back(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) { - + if (electrostaticAtomData.is_Fluctuating) { + ret = FQtypes.insert( atid ); + if (ret.second == false) { + sprintf( painCave.errMsg, + "Electrostatic already had a previous fluctuating charge entry with ident %d\n", + atid ); + painCave.severity = OPENMD_INFO; + painCave.isFatal = 0; + simError(); + } + FQtids[atid] = fqtid; + Jij[fqtid].resize(nFlucq_); + + // Now, iterate over all known fluctuating and add to the coulomb integral map: + + std::set::iterator it; + for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { + int etid2 = Etids[ (*it) ]; + int fqtid2 = FQtids[ (*it) ]; + ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ]; 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; @@ -709,16 +737,11 @@ namespace OpenMD { CubicSpline* J = new CubicSpline(); J->addPoints(rvals, Jvals); - - pair key1, key2; - key1 = make_pair(atomType, atype2); - key2 = make_pair(atype2, atomType); - - Jij[key1] = J; - Jij[key2] = J; - } - } - + Jij[fqtid][fqtid2] = J; + Jij[fqtid2].resize( nFlucq_ ); + Jij[fqtid2][fqtid] = J; + } + } return; } @@ -744,67 +767,46 @@ namespace OpenMD { void Electrostatic::calcForce(InteractionData &idat) { - RealType C_a, C_b; // Charges - Vector3d D_a, D_b; // Dipoles (space-fixed) - Mat3x3d Q_a, Q_b; // Quadrupoles (space-fixed) + if (!initialized_) initialize(); + + data1 = ElectrostaticMap[Etids[idat.atid1]]; + data2 = ElectrostaticMap[Etids[idat.atid2]]; - RealType ri; // Distance utility scalar - RealType rdDa, rdDb; // Dipole utility scalars - Vector3d rxDa, rxDb; // Dipole utility vectors - RealType rdQar, rdQbr, trQa, trQb; // Quadrupole utility scalars - Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr; // Quadrupole utility vectors - RealType pref; - - RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars - RealType rQaQbr; - Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors - Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; - Mat3x3d QaQb; // Cross-interaction matrices - - RealType U(0.0); // Potential - Vector3d F(0.0); // Force - Vector3d Ta(0.0); // Torque on site a - Vector3d Tb(0.0); // Torque on site b - Vector3d Ea(0.0); // Electric field at site a - Vector3d Eb(0.0); // Electric field at site b - RealType dUdCa(0.0); // fluctuating charge force at site a - RealType dUdCb(0.0); // fluctuating charge force at site a + U = 0.0; // Potential + F.zero(); // Force + Ta.zero(); // Torque on site a + Tb.zero(); // Torque on site b + Ea.zero(); // Electric field at site a + Eb.zero(); // Electric field at site b + dUdCa = 0.0; // fluctuating charge force at site a + dUdCb = 0.0; // fluctuating charge force at site a // Indirect interactions mediated by the reaction field. - RealType indirect_Pot(0.0); // Potential - Vector3d indirect_F(0.0); // Force - Vector3d indirect_Ta(0.0); // Torque on site a - Vector3d indirect_Tb(0.0); // Torque on site b + indirect_Pot = 0.0; // Potential + indirect_F.zero(); // Force + indirect_Ta.zero(); // Torque on site a + indirect_Tb.zero(); // Torque on site b // Excluded potential that is still computed for fluctuating charges - RealType excluded_Pot(0.0); + excluded_Pot= 0.0; - RealType rfContrib, coulInt; - - // spline for coulomb integral - CubicSpline* J; - if (!initialized_) initialize(); - - ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; - ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second]; - // some variables we'll need independent of electrostatic type: ri = 1.0 / *(idat.rij); - Vector3d rhat = *(idat.d) * ri; + rhat = *(idat.d) * ri; // logicals - bool a_is_Charge = data1.is_Charge; - bool a_is_Dipole = data1.is_Dipole; - bool a_is_Quadrupole = data1.is_Quadrupole; - bool a_is_Fluctuating = data1.is_Fluctuating; + a_is_Charge = data1.is_Charge; + a_is_Dipole = data1.is_Dipole; + a_is_Quadrupole = data1.is_Quadrupole; + a_is_Fluctuating = data1.is_Fluctuating; - bool b_is_Charge = data2.is_Charge; - bool b_is_Dipole = data2.is_Dipole; - bool b_is_Quadrupole = data2.is_Quadrupole; - bool b_is_Fluctuating = data2.is_Fluctuating; + b_is_Charge = data2.is_Charge; + b_is_Dipole = data2.is_Dipole; + b_is_Quadrupole = data2.is_Quadrupole; + b_is_Fluctuating = data2.is_Fluctuating; // Obtain all of the required radial function values from the // spline structures: @@ -911,7 +913,7 @@ namespace OpenMD { } if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { - J = Jij[idat.atypes]; + J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; } if (a_is_Charge) { @@ -1149,7 +1151,7 @@ namespace OpenMD { if (!initialized_) initialize(); - ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; + ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]]; // logicals bool i_is_Charge = data.is_Charge;