| 53 |
|
#include "io/Globals.hpp" |
| 54 |
|
#include "nonbonded/SlaterIntegrals.hpp" |
| 55 |
|
#include "utils/PhysicalConstants.hpp" |
| 56 |
+ |
#include "math/erfc.hpp" |
| 57 |
|
|
| 57 |
– |
|
| 58 |
|
namespace OpenMD { |
| 59 |
|
|
| 60 |
|
Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), |
| 361 |
|
vector<RealType> rvals; |
| 362 |
|
vector<RealType> J1vals; |
| 363 |
|
vector<RealType> J2vals; |
| 364 |
< |
for (int i = 0; i < np_; i++) { |
| 364 |
> |
// don't start at i = 0, as rval = 0 is undefined for the slater overlap integrals. |
| 365 |
> |
for (int i = 1; i < np_+1; i++) { |
| 366 |
|
rval = RealType(i) * dr; |
| 367 |
|
rvals.push_back(rval); |
| 368 |
< |
J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); |
| 368 |
> |
J1vals.push_back(sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal ); |
| 369 |
|
// may not be necessary if Slater coulomb integral is symmetric |
| 370 |
< |
J2vals.push_back(eaData2.hardness * sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); |
| 370 |
> |
J2vals.push_back(sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromToBohr ) * PhysicalConstants::hartreeToKcal ); |
| 371 |
|
} |
| 372 |
|
|
| 373 |
|
CubicSpline* J1 = new CubicSpline(); |
| 579 |
|
if (j_is_Charge) { |
| 580 |
|
if (screeningMethod_ == DAMPED) { |
| 581 |
|
// assemble the damping variables |
| 582 |
< |
//res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); |
| 583 |
< |
//erfcVal = res.first; |
| 584 |
< |
//derfcVal = res.second; |
| 582 |
> |
res = erfcSpline_->getValueAndDerivativeAt( *(idat.rij) ); |
| 583 |
> |
erfcVal = res.first; |
| 584 |
> |
derfcVal = res.second; |
| 585 |
|
|
| 586 |
< |
erfcVal = erfc(dampingAlpha_ * *(idat.rij)); |
| 587 |
< |
derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); |
| 586 |
> |
//erfcVal = erfc(dampingAlpha_ * *(idat.rij)); |
| 587 |
> |
//derfcVal = - alphaPi_ * exp(-alpha2_ * *(idat.r2)); |
| 588 |
|
|
| 589 |
|
c1 = erfcVal * riji; |
| 590 |
|
c2 = (-derfcVal + c1) * riji; |
| 1050 |
|
// indirect reaction field terms. |
| 1051 |
|
|
| 1052 |
|
*(idat.vpair) += indirect_vpair; |
| 1053 |
+ |
|
| 1054 |
+ |
(*(idat.excludedPot))[ELECTROSTATIC_FAMILY] += (*(idat.sw) * vterm + |
| 1055 |
+ |
vFluc1 ) * q_i * q_j; |
| 1056 |
|
(*(idat.pot))[ELECTROSTATIC_FAMILY] += indirect_Pot; |
| 1057 |
|
*(idat.f1) += indirect_dVdr; |
| 1058 |
|
|
| 1081 |
|
chg1 += *(sdat.flucQ); |
| 1082 |
|
// dVdFQ is really a force, so this is negative the derivative |
| 1083 |
|
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1084 |
+ |
(*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) * |
| 1085 |
+ |
(*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity); |
| 1086 |
|
} |
| 1087 |
|
|
| 1088 |
|
if (summationMethod_ == esm_REACTION_FIELD) { |
