| 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) { |
| 637 |
|
} |
| 638 |
|
|
| 639 |
|
void Electrostatic::addType(AtomType* atomType){ |
| 640 |
< |
|
| 640 |
> |
|
| 641 |
|
ElectrostaticAtomData electrostaticAtomData; |
| 642 |
|
electrostaticAtomData.is_Charge = false; |
| 643 |
|
electrostaticAtomData.is_Dipole = false; |
| 673 |
|
electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); |
| 674 |
|
} |
| 675 |
|
|
| 676 |
< |
pair<map<int,AtomType*>::iterator,bool> ret; |
| 677 |
< |
ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(), |
| 678 |
< |
atomType) ); |
| 676 |
> |
int atid = atomType->getIdent(); |
| 677 |
> |
int etid = Etypes.size(); |
| 678 |
> |
int fqtid = FQtypes.size(); |
| 679 |
> |
|
| 680 |
> |
pair<set<int>::iterator,bool> ret; |
| 681 |
> |
ret = Etypes.insert( atid ); |
| 682 |
|
if (ret.second == false) { |
| 683 |
|
sprintf( painCave.errMsg, |
| 684 |
|
"Electrostatic already had a previous entry with ident %d\n", |
| 685 |
< |
atomType->getIdent() ); |
| 685 |
> |
atid); |
| 686 |
|
painCave.severity = OPENMD_INFO; |
| 687 |
|
painCave.isFatal = 0; |
| 688 |
|
simError(); |
| 689 |
|
} |
| 690 |
|
|
| 691 |
< |
ElectrostaticMap[atomType] = electrostaticAtomData; |
| 691 |
> |
Etids[ atid ] = etid; |
| 692 |
> |
ElectrostaticMap.push_back(electrostaticAtomData); |
| 693 |
|
|
| 694 |
< |
// Now, iterate over all known types and add to the mixing map: |
| 695 |
< |
|
| 696 |
< |
map<AtomType*, ElectrostaticAtomData>::iterator it; |
| 697 |
< |
for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) { |
| 698 |
< |
AtomType* atype2 = (*it).first; |
| 699 |
< |
ElectrostaticAtomData eaData2 = (*it).second; |
| 700 |
< |
if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) { |
| 701 |
< |
|
| 694 |
> |
if (electrostaticAtomData.is_Fluctuating) { |
| 695 |
> |
ret = FQtypes.insert( atid ); |
| 696 |
> |
if (ret.second == false) { |
| 697 |
> |
sprintf( painCave.errMsg, |
| 698 |
> |
"Electrostatic already had a previous fluctuating charge entry with ident %d\n", |
| 699 |
> |
atid ); |
| 700 |
> |
painCave.severity = OPENMD_INFO; |
| 701 |
> |
painCave.isFatal = 0; |
| 702 |
> |
simError(); |
| 703 |
> |
} |
| 704 |
> |
FQtids[atid] = fqtid; |
| 705 |
> |
Jij[fqtid].resize(nFlucq_); |
| 706 |
> |
|
| 707 |
> |
// Now, iterate over all known fluctuating and add to the coulomb integral map: |
| 708 |
> |
|
| 709 |
> |
std::set<int>::iterator it; |
| 710 |
> |
for( it = FQtypes.begin(); it != FQtypes.end(); ++it) { |
| 711 |
> |
int etid2 = Etids[ (*it) ]; |
| 712 |
> |
int fqtid2 = FQtids[ (*it) ]; |
| 713 |
> |
ElectrostaticAtomData eaData2 = ElectrostaticMap[ etid2 ]; |
| 714 |
|
RealType a = electrostaticAtomData.slaterZeta; |
| 715 |
|
RealType b = eaData2.slaterZeta; |
| 716 |
|
int m = electrostaticAtomData.slaterN; |
| 717 |
|
int n = eaData2.slaterN; |
| 718 |
< |
|
| 718 |
> |
|
| 719 |
|
// Create the spline of the coulombic integral for s-type |
| 720 |
|
// Slater orbitals. Add a 2 angstrom safety window to deal |
| 721 |
|
// with cutoffGroups that have charged atoms longer than the |
| 722 |
|
// cutoffRadius away from each other. |
| 723 |
< |
|
| 723 |
> |
|
| 724 |
|
RealType rval; |
| 725 |
|
RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1); |
| 726 |
|
vector<RealType> rvals; |
| 737 |
|
|
| 738 |
|
CubicSpline* J = new CubicSpline(); |
| 739 |
|
J->addPoints(rvals, Jvals); |
| 740 |
< |
|
| 741 |
< |
pair<AtomType*, AtomType*> key1, key2; |
| 742 |
< |
key1 = make_pair(atomType, atype2); |
| 743 |
< |
key2 = make_pair(atype2, atomType); |
| 744 |
< |
|
| 717 |
< |
Jij[key1] = J; |
| 718 |
< |
Jij[key2] = J; |
| 719 |
< |
} |
| 720 |
< |
} |
| 721 |
< |
|
| 740 |
> |
Jij[fqtid][fqtid2] = J; |
| 741 |
> |
Jij[fqtid2].resize( nFlucq_ ); |
| 742 |
> |
Jij[fqtid2][fqtid] = J; |
| 743 |
> |
} |
| 744 |
> |
} |
| 745 |
|
return; |
| 746 |
|
} |
| 747 |
|
|
| 769 |
|
|
| 770 |
|
if (!initialized_) initialize(); |
| 771 |
|
|
| 772 |
< |
data1 = ElectrostaticMap[idat.atypes.first]; |
| 773 |
< |
data2 = ElectrostaticMap[idat.atypes.second]; |
| 772 |
> |
data1 = ElectrostaticMap[Etids[idat.atid1]]; |
| 773 |
> |
data2 = ElectrostaticMap[Etids[idat.atid2]]; |
| 774 |
|
|
| 775 |
|
U = 0.0; // Potential |
| 776 |
|
F.zero(); // Force |
| 913 |
|
} |
| 914 |
|
|
| 915 |
|
if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) { |
| 916 |
< |
J = Jij[idat.atypes]; |
| 916 |
> |
J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]]; |
| 917 |
|
} |
| 918 |
|
|
| 919 |
|
if (a_is_Charge) { |
| 1151 |
|
|
| 1152 |
|
if (!initialized_) initialize(); |
| 1153 |
|
|
| 1154 |
< |
ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; |
| 1154 |
> |
ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]]; |
| 1155 |
|
|
| 1156 |
|
// logicals |
| 1157 |
|
bool i_is_Charge = data.is_Charge; |
