--- trunk/src/nonbonded/Electrostatic.cpp 2013/06/16 15:15:42 1879 +++ trunk/src/nonbonded/Electrostatic.cpp 2013/07/19 18:18:27 1907 @@ -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) { @@ -250,7 +262,13 @@ namespace OpenMD { b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; - selfMult_ = b0c + a2 * invArootPi; + //selfMult1_ = - 2.0 * a2 * invArootPi; + //selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0; + //selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0; + // Half the Smith self piece: + selfMult1_ = - a2 * invArootPi; + selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0; + selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0; } else { a2 = 0.0; b0c = 1.0 / r; @@ -259,7 +277,9 @@ namespace OpenMD { b3c = (5.0 * b2c) / r2; b4c = (7.0 * b3c) / r2; b5c = (9.0 * b4c) / r2; - selfMult_ = b0c; + selfMult1_ = 0.0; + selfMult2_ = 0.0; + selfMult4_ = 0.0; } // higher derivatives of B_0 at the cutoff radius: @@ -267,8 +287,11 @@ namespace OpenMD { db0c_2 = -b1c + r2 * b2c; db0c_3 = 3.0*r*b2c - r2*r*b3c; db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; - db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; + db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; + selfMult1_ -= b0c; + selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; + selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; // working variables for the splines: RealType ri, ri2; @@ -398,7 +421,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 +478,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 +538,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 +648,7 @@ namespace OpenMD { } void Electrostatic::addType(AtomType* atomType){ - + ElectrostaticAtomData electrostaticAtomData; electrostaticAtomData.is_Charge = false; electrostaticAtomData.is_Dipole = false; @@ -661,38 +684,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 +748,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 +778,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 +924,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) { @@ -1102,7 +1115,6 @@ namespace OpenMD { Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; - // cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; } } @@ -1149,15 +1161,20 @@ namespace OpenMD { if (!initialized_) initialize(); - ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; + ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]]; // logicals bool i_is_Charge = data.is_Charge; bool i_is_Dipole = data.is_Dipole; + bool i_is_Quadrupole = data.is_Quadrupole; bool i_is_Fluctuating = data.is_Fluctuating; RealType C_a = data.fixedCharge; - RealType self, preVal, DadDa; - + RealType self(0.0), preVal, DdD, trQ, trQQ; + + if (i_is_Dipole) { + DdD = data.dipole.lengthSquare(); + } + if (i_is_Fluctuating) { C_a += *(sdat.flucQ); // dVdFQ is really a force, so this is negative the derivative @@ -1178,18 +1195,26 @@ namespace OpenMD { } if (i_is_Dipole) { - DadDa = data.dipole.lengthSquare(); - (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa; + (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD; } break; case esm_SHIFTED_FORCE: case esm_SHIFTED_POTENTIAL: - if (i_is_Charge) { - self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; - (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; + case esm_TAYLOR_SHIFTED: + if (i_is_Charge) + self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); + if (i_is_Dipole) + self += selfMult2_ * pre22_ * DdD; + if (i_is_Quadrupole) { + trQ = data.quadrupole.trace(); + trQQ = (data.quadrupole * data.quadrupole).trace(); + self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ); + if (i_is_Charge) + self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ; } + (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; break; default: break;