| 213 |
|
haveDampingAlpha_ = true; |
| 214 |
|
} |
| 215 |
|
|
| 216 |
< |
// find all of the Electrostatic atom Types: |
| 217 |
< |
ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes(); |
| 218 |
< |
ForceField::AtomTypeContainer::MapTypeIterator i; |
| 219 |
< |
AtomType* at; |
| 216 |
> |
Etypes.clear(); |
| 217 |
> |
Etids.clear(); |
| 218 |
> |
FQtypes.clear(); |
| 219 |
> |
FQtids.clear(); |
| 220 |
> |
ElectrostaticMap.clear(); |
| 221 |
> |
Jij.clear(); |
| 222 |
> |
nElectro_ = 0; |
| 223 |
> |
nFlucq_ = 0; |
| 224 |
> |
|
| 225 |
> |
Etids.resize( forceField_->getNAtomType(), -1); |
| 226 |
> |
FQtids.resize( forceField_->getNAtomType(), -1); |
| 227 |
> |
|
| 228 |
> |
set<AtomType*>::iterator at; |
| 229 |
> |
for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { |
| 230 |
> |
if ((*at)->isElectrostatic()) nElectro_++; |
| 231 |
> |
if ((*at)->isFluctuatingCharge()) nFlucq_++; |
| 232 |
> |
} |
| 233 |
|
|
| 234 |
< |
for (at = atomTypes->beginType(i); at != NULL; |
| 235 |
< |
at = atomTypes->nextType(i)) { |
| 236 |
< |
|
| 237 |
< |
if (at->isElectrostatic()) |
| 225 |
< |
addType(at); |
| 234 |
> |
Jij.resize(nFlucq_); |
| 235 |
> |
|
| 236 |
> |
for (at = simTypes_.begin(); at != simTypes_.end(); ++at) { |
| 237 |
> |
if ((*at)->isElectrostatic()) addType(*at); |
| 238 |
|
} |
| 239 |
|
|
| 240 |
|
if (summationMethod_ == esm_REACTION_FIELD) { |
| 262 |
|
b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2; |
| 263 |
|
b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2; |
| 264 |
|
b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2; |
| 265 |
< |
selfMult_ = b0c + a2 * invArootPi; |
| 265 |
> |
//selfMult1_ = - 2.0 * a2 * invArootPi; |
| 266 |
> |
//selfMult2_ = - 4.0 * a2 * a2 * invArootPi / 3.0; |
| 267 |
> |
//selfMult4_ = - 8.0 * a2 * a2 * a2 * invArootPi / 5.0; |
| 268 |
> |
// Half the Smith self piece: |
| 269 |
> |
selfMult1_ = - a2 * invArootPi; |
| 270 |
> |
selfMult2_ = - 2.0 * a2 * a2 * invArootPi / 3.0; |
| 271 |
> |
selfMult4_ = - 4.0 * a2 * a2 * a2 * invArootPi / 5.0; |
| 272 |
|
} else { |
| 273 |
|
a2 = 0.0; |
| 274 |
|
b0c = 1.0 / r; |
| 277 |
|
b3c = (5.0 * b2c) / r2; |
| 278 |
|
b4c = (7.0 * b3c) / r2; |
| 279 |
|
b5c = (9.0 * b4c) / r2; |
| 280 |
< |
selfMult_ = b0c; |
| 280 |
> |
selfMult1_ = 0.0; |
| 281 |
> |
selfMult2_ = 0.0; |
| 282 |
> |
selfMult4_ = 0.0; |
| 283 |
|
} |
| 284 |
|
|
| 285 |
|
// higher derivatives of B_0 at the cutoff radius: |
| 287 |
|
db0c_2 = -b1c + r2 * b2c; |
| 288 |
|
db0c_3 = 3.0*r*b2c - r2*r*b3c; |
| 289 |
|
db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c; |
| 290 |
< |
db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; |
| 290 |
> |
db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c; |
| 291 |
|
|
| 292 |
+ |
selfMult1_ -= b0c; |
| 293 |
+ |
selfMult2_ += (db0c_2 + 2.0*db0c_1*ric) / 3.0; |
| 294 |
+ |
selfMult4_ -= (db0c_4 + 4.0*db0c_3*ric) / 15.0; |
| 295 |
|
|
| 296 |
|
// working variables for the splines: |
| 297 |
|
RealType ri, ri2; |
| 648 |
|
} |
| 649 |
|
|
| 650 |
|
void Electrostatic::addType(AtomType* atomType){ |
| 651 |
< |
|
| 651 |
> |
|
| 652 |
|
ElectrostaticAtomData electrostaticAtomData; |
| 653 |
|
electrostaticAtomData.is_Charge = false; |
| 654 |
|
electrostaticAtomData.is_Dipole = false; |
| 684 |
|
electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); |
| 685 |
|
} |
| 686 |
|
|
| 687 |
< |
pair<map<int,AtomType*>::iterator,bool> ret; |
| 688 |
< |
ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(), |
| 689 |
< |
atomType) ); |
| 687 |
> |
int atid = atomType->getIdent(); |
| 688 |
> |
int etid = Etypes.size(); |
| 689 |
> |
int fqtid = FQtypes.size(); |
| 690 |
> |
|
| 691 |
> |
pair<set<int>::iterator,bool> ret; |
| 692 |
> |
ret = Etypes.insert( atid ); |
| 693 |
|
if (ret.second == false) { |
| 694 |
|
sprintf( painCave.errMsg, |
| 695 |
|
"Electrostatic already had a previous entry with ident %d\n", |
| 696 |
< |
atomType->getIdent() ); |
| 696 |
> |
atid); |
| 697 |
|
painCave.severity = OPENMD_INFO; |
| 698 |
|
painCave.isFatal = 0; |
| 699 |
|
simError(); |
| 700 |
|
} |
| 701 |
|
|
| 702 |
< |
ElectrostaticMap[atomType] = electrostaticAtomData; |
| 702 |
> |
Etids[ atid ] = etid; |
| 703 |
> |
ElectrostaticMap.push_back(electrostaticAtomData); |
| 704 |
|
|
| 705 |
< |
// Now, iterate over all known types and add to the mixing map: |
| 706 |
< |
|
| 707 |
< |
map<AtomType*, ElectrostaticAtomData>::iterator it; |
| 708 |
< |
for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) { |
| 709 |
< |
AtomType* atype2 = (*it).first; |
| 710 |
< |
ElectrostaticAtomData eaData2 = (*it).second; |
| 711 |
< |
if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) { |
| 712 |
< |
|
| 705 |
> |
if (electrostaticAtomData.is_Fluctuating) { |
| 706 |
> |
ret = FQtypes.insert( atid ); |
| 707 |
> |
if (ret.second == false) { |
| 708 |
> |
sprintf( painCave.errMsg, |
| 709 |
> |
"Electrostatic already had a previous fluctuating charge entry with ident %d\n", |
| 710 |
> |
atid ); |
| 711 |
> |
painCave.severity = OPENMD_INFO; |
| 712 |
> |
painCave.isFatal = 0; |
| 713 |
> |
simError(); |
| 714 |
> |
} |
| 715 |
> |
FQtids[atid] = fqtid; |
| 716 |
> |
Jij[fqtid].resize(nFlucq_); |
| 717 |
> |
|
| 718 |
> |
// Now, iterate over all known fluctuating and add to the coulomb integral map: |
| 719 |
> |
|
| 720 |
> |
std::set<int>::iterator it; |
| 721 |
> |
for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { |
| 722 |
> |
int etid2 = Etids[ (*it) ]; |
| 723 |
> |
int fqtid2 = FQtids[ (*it) ]; |
| 724 |
> |
ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ]; |
| 725 |
|
RealType a = electrostaticAtomData.slaterZeta; |
| 726 |
|
RealType b = eaData2.slaterZeta; |
| 727 |
|
int m = electrostaticAtomData.slaterN; |
| 728 |
|
int n = eaData2.slaterN; |
| 729 |
< |
|
| 729 |
> |
|
| 730 |
|
// Create the spline of the coulombic integral for s-type |
| 731 |
|
// Slater orbitals. Add a 2 angstrom safety window to deal |
| 732 |
|
// with cutoffGroups that have charged atoms longer than the |
| 733 |
|
// cutoffRadius away from each other. |
| 734 |
< |
|
| 734 |
> |
|
| 735 |
|
RealType rval; |
| 736 |
|
RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1); |
| 737 |
|
vector<RealType> rvals; |
| 748 |
|
|
| 749 |
|
CubicSpline* J = new CubicSpline(); |
| 750 |
|
J->addPoints(rvals, Jvals); |
| 751 |
< |
|
| 752 |
< |
pair<AtomType*, AtomType*> key1, key2; |
| 753 |
< |
key1 = make_pair(atomType, atype2); |
| 754 |
< |
key2 = make_pair(atype2, atomType); |
| 755 |
< |
|
| 717 |
< |
Jij[key1] = J; |
| 718 |
< |
Jij[key2] = J; |
| 719 |
< |
} |
| 720 |
< |
} |
| 721 |
< |
|
| 751 |
> |
Jij[fqtid][fqtid2] = J; |
| 752 |
> |
Jij[fqtid2].resize( nFlucq_ ); |
| 753 |
> |
Jij[fqtid2][fqtid] = J; |
| 754 |
> |
} |
| 755 |
> |
} |
| 756 |
|
return; |
| 757 |
|
} |
| 758 |
|
|
| 780 |
|
|
| 781 |
|
if (!initialized_) initialize(); |
| 782 |
|
|
| 783 |
< |
data1 = ElectrostaticMap[idat.atypes.first]; |
| 784 |
< |
data2 = ElectrostaticMap[idat.atypes.second]; |
| 783 |
> |
data1 = ElectrostaticMap[Etids[idat.atid1]]; |
| 784 |
> |
data2 = ElectrostaticMap[Etids[idat.atid2]]; |
| 785 |
|
|
| 786 |
|
U = 0.0; // Potential |
| 787 |
|
F.zero(); // Force |
| 924 |
|
} |
| 925 |
|
|
| 926 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 927 |
< |
J = Jij[idat.atypes]; |
| 927 |
> |
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 928 |
|
} |
| 929 |
|
|
| 930 |
|
if (a_is_Charge) { |
| 1115 |
|
|
| 1116 |
|
Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43; |
| 1117 |
|
|
| 1084 |
– |
// cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n"; |
| 1118 |
|
} |
| 1119 |
|
} |
| 1120 |
|
|
| 1161 |
|
|
| 1162 |
|
if (!initialized_) initialize(); |
| 1163 |
|
|
| 1164 |
< |
ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; |
| 1164 |
> |
ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]]; |
| 1165 |
|
|
| 1166 |
|
// logicals |
| 1167 |
|
bool i_is_Charge = data.is_Charge; |
| 1168 |
|
bool i_is_Dipole = data.is_Dipole; |
| 1169 |
+ |
bool i_is_Quadrupole = data.is_Quadrupole; |
| 1170 |
|
bool i_is_Fluctuating = data.is_Fluctuating; |
| 1171 |
|
RealType C_a = data.fixedCharge; |
| 1172 |
< |
RealType self, preVal, DadDa; |
| 1173 |
< |
|
| 1172 |
> |
RealType self(0.0), preVal, DdD, trQ, trQQ; |
| 1173 |
> |
|
| 1174 |
> |
if (i_is_Dipole) { |
| 1175 |
> |
DdD = data.dipole.lengthSquare(); |
| 1176 |
> |
} |
| 1177 |
> |
|
| 1178 |
|
if (i_is_Fluctuating) { |
| 1179 |
|
C_a += *(sdat.flucQ); |
| 1180 |
|
// dVdFQ is really a force, so this is negative the derivative |
| 1195 |
|
} |
| 1196 |
|
|
| 1197 |
|
if (i_is_Dipole) { |
| 1198 |
< |
DadDa = data.dipole.lengthSquare(); |
| 1161 |
< |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa; |
| 1198 |
> |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DdD; |
| 1199 |
|
} |
| 1200 |
|
|
| 1201 |
|
break; |
| 1202 |
|
|
| 1203 |
|
case esm_SHIFTED_FORCE: |
| 1204 |
|
case esm_SHIFTED_POTENTIAL: |
| 1205 |
< |
if (i_is_Charge) { |
| 1206 |
< |
self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_; |
| 1207 |
< |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; |
| 1205 |
> |
case esm_TAYLOR_SHIFTED: |
| 1206 |
> |
if (i_is_Charge) |
| 1207 |
> |
self += selfMult1_ * pre11_ * C_a * (C_a + *(sdat.skippedCharge)); |
| 1208 |
> |
if (i_is_Dipole) |
| 1209 |
> |
self += selfMult2_ * pre22_ * DdD; |
| 1210 |
> |
if (i_is_Quadrupole) { |
| 1211 |
> |
trQ = data.quadrupole.trace(); |
| 1212 |
> |
trQQ = (data.quadrupole * data.quadrupole).trace(); |
| 1213 |
> |
self += selfMult4_ * pre44_ * (2.0*trQQ + trQ*trQ); |
| 1214 |
> |
if (i_is_Charge) |
| 1215 |
> |
self -= selfMult2_ * pre14_ * 2.0 * C_a * trQ; |
| 1216 |
|
} |
| 1217 |
+ |
(*(sdat.pot))[ELECTROSTATIC_FAMILY] += self; |
| 1218 |
|
break; |
| 1219 |
|
default: |
| 1220 |
|
break; |
