| 50 |
|
#include "types/FixedChargeAdapter.hpp" |
| 51 |
|
#include "types/MultipoleAdapter.hpp" |
| 52 |
|
#include "io/Globals.hpp" |
| 53 |
+ |
#include "nonbonded/SlaterIntegrals.hpp" |
| 54 |
+ |
#include "utils/PhysicalConstants.hpp" |
| 55 |
|
|
| 56 |
+ |
|
| 57 |
|
namespace OpenMD { |
| 58 |
|
|
| 59 |
|
Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), |
| 312 |
|
} |
| 313 |
|
} |
| 314 |
|
|
| 315 |
+ |
FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); |
| 316 |
+ |
|
| 317 |
+ |
if (fqa.isFluctuatingCharge()) { |
| 318 |
+ |
electrostaticAtomData.is_FluctuatingCharge = true; |
| 319 |
+ |
electrostaticAtomData.electronegativity = fca.getElectronegativity(); |
| 320 |
+ |
electrostaticAtomData.hardness = fca.getHardness(); |
| 321 |
+ |
electrostaticAtomData.slaterN = fca.getSlaterN(); |
| 322 |
+ |
electrostaticAtomData.slaterZeta = fca.getSlaterZeta(); |
| 323 |
+ |
} |
| 324 |
|
|
| 325 |
|
pair<map<int,AtomType*>::iterator,bool> ret; |
| 326 |
|
ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(), |
| 334 |
|
simError(); |
| 335 |
|
} |
| 336 |
|
|
| 337 |
< |
ElectrostaticMap[atomType] = electrostaticAtomData; |
| 337 |
> |
ElectrostaticMap[atomType] = electrostaticAtomData; |
| 338 |
> |
|
| 339 |
> |
// Now, iterate over all known types and add to the mixing map: |
| 340 |
> |
|
| 341 |
> |
map<AtomType*, ElectrostaticAtomData>::iterator it; |
| 342 |
> |
for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) { |
| 343 |
> |
AtomType* atype2 = (*it).first; |
| 344 |
> |
|
| 345 |
> |
if ((*it).is_FluctuatingCharge && electrostaticAtomData.is_FluctuatingCharge) { |
| 346 |
> |
|
| 347 |
> |
RealType a = electrostaticAtomData.slaterZeta; |
| 348 |
> |
RealType b = (*it).slaterZeta; |
| 349 |
> |
int m = electrostaticAtomData.slaterN; |
| 350 |
> |
int n = (*it).slaterN; |
| 351 |
> |
|
| 352 |
> |
// Create the spline of the coulombic integral for s-type |
| 353 |
> |
// Slater orbitals. Add a 2 angstrom safety window to deal |
| 354 |
> |
// with cutoffGroups that have charged atoms longer than the |
| 355 |
> |
// cutoffRadius away from each other. |
| 356 |
> |
|
| 357 |
> |
RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1); |
| 358 |
> |
vector<RealType> rvals; |
| 359 |
> |
vector<RealType> J1vals; |
| 360 |
> |
vector<RealType> J2vals; |
| 361 |
> |
for (int i = 0; i < np_; i++) { |
| 362 |
> |
rval = RealType(i) * dr; |
| 363 |
> |
rvals.push_back(rval); |
| 364 |
> |
J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); |
| 365 |
> |
J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); |
| 366 |
> |
} |
| 367 |
> |
|
| 368 |
> |
CubicSpline J1 = new CubicSpline(); |
| 369 |
> |
J1->addPoints(rvals, J1vals); |
| 370 |
> |
CubicSpline J2 = new CubicSpline(); |
| 371 |
> |
J2->addPoints(rvals, J2vals); |
| 372 |
> |
|
| 373 |
> |
pair<AtomType*, AtomType*> key1, key2; |
| 374 |
> |
key1 = make_pair(atomType, atype2); |
| 375 |
> |
key2 = make_pair(atype2, atomType); |
| 376 |
> |
|
| 377 |
> |
Jij[key1] = J1; |
| 378 |
> |
Jij[key2] = J2; |
| 379 |
> |
} |
| 380 |
> |
} |
| 381 |
> |
|
| 382 |
|
return; |
| 383 |
|
} |
| 384 |
|
|
