| 48 |
|
#include "utils/simError.h" |
| 49 |
|
#include "types/NonBondedInteractionType.hpp" |
| 50 |
|
#include "types/FixedChargeAdapter.hpp" |
| 51 |
+ |
#include "types/FluctuatingChargeAdapter.hpp" |
| 52 |
|
#include "types/MultipoleAdapter.hpp" |
| 53 |
|
#include "io/Globals.hpp" |
| 54 |
+ |
#include "nonbonded/SlaterIntegrals.hpp" |
| 55 |
+ |
#include "utils/PhysicalConstants.hpp" |
| 56 |
|
|
| 57 |
+ |
|
| 58 |
|
namespace OpenMD { |
| 59 |
|
|
| 60 |
|
Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false), |
| 193 |
|
|
| 194 |
|
// throw warning |
| 195 |
|
sprintf( painCave.errMsg, |
| 196 |
< |
"Electrostatic::initialize: dampingAlpha was not specified in the input file.\n" |
| 197 |
< |
"\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", |
| 196 |
> |
"Electrostatic::initialize: dampingAlpha was not specified in the\n" |
| 197 |
> |
"\tinput file. A default value of %f (1/ang) will be used for the\n" |
| 198 |
> |
"\tcutoff of %f (ang).\n", |
| 199 |
|
dampingAlpha_, cutoffRadius_); |
| 200 |
|
painCave.severity = OPENMD_INFO; |
| 201 |
|
painCave.isFatal = 0; |
| 218 |
|
addType(at); |
| 219 |
|
} |
| 220 |
|
|
| 216 |
– |
|
| 221 |
|
cutoffRadius2_ = cutoffRadius_ * cutoffRadius_; |
| 222 |
|
rcuti_ = 1.0 / cutoffRadius_; |
| 223 |
|
rcuti2_ = rcuti_ * rcuti_; |
| 284 |
|
electrostaticAtomData.is_Dipole = false; |
| 285 |
|
electrostaticAtomData.is_SplitDipole = false; |
| 286 |
|
electrostaticAtomData.is_Quadrupole = false; |
| 287 |
+ |
electrostaticAtomData.is_Fluctuating = false; |
| 288 |
|
|
| 289 |
|
FixedChargeAdapter fca = FixedChargeAdapter(atomType); |
| 290 |
|
|
| 291 |
|
if (fca.isFixedCharge()) { |
| 292 |
|
electrostaticAtomData.is_Charge = true; |
| 293 |
< |
electrostaticAtomData.charge = fca.getCharge(); |
| 293 |
> |
electrostaticAtomData.fixedCharge = fca.getCharge(); |
| 294 |
|
} |
| 295 |
|
|
| 296 |
|
MultipoleAdapter ma = MultipoleAdapter(atomType); |
| 314 |
|
} |
| 315 |
|
} |
| 316 |
|
|
| 317 |
+ |
FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); |
| 318 |
|
|
| 319 |
+ |
if (fqa.isFluctuatingCharge()) { |
| 320 |
+ |
electrostaticAtomData.is_Fluctuating = true; |
| 321 |
+ |
electrostaticAtomData.electronegativity = fqa.getElectronegativity(); |
| 322 |
+ |
electrostaticAtomData.hardness = fqa.getHardness(); |
| 323 |
+ |
electrostaticAtomData.slaterN = fqa.getSlaterN(); |
| 324 |
+ |
electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); |
| 325 |
+ |
} |
| 326 |
+ |
|
| 327 |
|
pair<map<int,AtomType*>::iterator,bool> ret; |
| 328 |
|
ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(), |
| 329 |
|
atomType) ); |
| 336 |
|
simError(); |
| 337 |
|
} |
| 338 |
|
|
| 339 |
< |
ElectrostaticMap[atomType] = electrostaticAtomData; |
| 339 |
> |
ElectrostaticMap[atomType] = electrostaticAtomData; |
| 340 |
> |
|
| 341 |
> |
// Now, iterate over all known types and add to the mixing map: |
| 342 |
> |
|
| 343 |
> |
map<AtomType*, ElectrostaticAtomData>::iterator it; |
| 344 |
> |
for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) { |
| 345 |
> |
AtomType* atype2 = (*it).first; |
| 346 |
> |
ElectrostaticAtomData eaData2 = (*it).second; |
| 347 |
> |
if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) { |
| 348 |
> |
|
| 349 |
> |
RealType a = electrostaticAtomData.slaterZeta; |
| 350 |
> |
RealType b = eaData2.slaterZeta; |
| 351 |
> |
int m = electrostaticAtomData.slaterN; |
| 352 |
> |
int n = eaData2.slaterN; |
| 353 |
> |
|
| 354 |
> |
// Create the spline of the coulombic integral for s-type |
| 355 |
> |
// Slater orbitals. Add a 2 angstrom safety window to deal |
| 356 |
> |
// with cutoffGroups that have charged atoms longer than the |
| 357 |
> |
// cutoffRadius away from each other. |
| 358 |
> |
|
| 359 |
> |
RealType rval; |
| 360 |
> |
RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1); |
| 361 |
> |
vector<RealType> rvals; |
| 362 |
> |
vector<RealType> J1vals; |
| 363 |
> |
vector<RealType> J2vals; |
| 364 |
> |
for (int i = 0; i < np_; i++) { |
| 365 |
> |
rval = RealType(i) * dr; |
| 366 |
> |
rvals.push_back(rval); |
| 367 |
> |
J1vals.push_back(electrostaticAtomData.hardness * sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); |
| 368 |
> |
// may not be necessary if Slater coulomb integral is symmetric |
| 369 |
> |
J2vals.push_back(eaData2.hardness * sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); |
| 370 |
> |
} |
| 371 |
> |
|
| 372 |
> |
CubicSpline* J1 = new CubicSpline(); |
| 373 |
> |
J1->addPoints(rvals, J1vals); |
| 374 |
> |
CubicSpline* J2 = new CubicSpline(); |
| 375 |
> |
J2->addPoints(rvals, J2vals); |
| 376 |
> |
|
| 377 |
> |
pair<AtomType*, AtomType*> key1, key2; |
| 378 |
> |
key1 = make_pair(atomType, atype2); |
| 379 |
> |
key2 = make_pair(atype2, atomType); |
| 380 |
> |
|
| 381 |
> |
Jij[key1] = J1; |
| 382 |
> |
Jij[key2] = J2; |
| 383 |
> |
} |
| 384 |
> |
} |
| 385 |
> |
|
| 386 |
|
return; |
| 387 |
|
} |
| 388 |
|
|
| 448 |
|
Vector3d indirect_dVdr(V3Zero); |
| 449 |
|
Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero); |
| 450 |
|
|
| 451 |
+ |
RealType coulInt, vFluc1(0.0), vFluc2(0.0); |
| 452 |
|
pair<RealType, RealType> res; |
| 453 |
|
|
| 454 |
+ |
// splines for coulomb integrals |
| 455 |
+ |
CubicSpline* J1; |
| 456 |
+ |
CubicSpline* J2; |
| 457 |
+ |
|
| 458 |
|
if (!initialized_) initialize(); |
| 459 |
|
|
| 460 |
|
ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; |
| 471 |
|
bool i_is_Dipole = data1.is_Dipole; |
| 472 |
|
bool i_is_SplitDipole = data1.is_SplitDipole; |
| 473 |
|
bool i_is_Quadrupole = data1.is_Quadrupole; |
| 474 |
+ |
bool i_is_Fluctuating = data1.is_Fluctuating; |
| 475 |
|
|
| 476 |
|
bool j_is_Charge = data2.is_Charge; |
| 477 |
|
bool j_is_Dipole = data2.is_Dipole; |
| 478 |
|
bool j_is_SplitDipole = data2.is_SplitDipole; |
| 479 |
|
bool j_is_Quadrupole = data2.is_Quadrupole; |
| 480 |
+ |
bool j_is_Fluctuating = data2.is_Fluctuating; |
| 481 |
|
|
| 482 |
|
if (i_is_Charge) { |
| 483 |
< |
q_i = data1.charge; |
| 483 |
> |
q_i = data1.fixedCharge; |
| 484 |
> |
|
| 485 |
> |
if (i_is_Fluctuating) { |
| 486 |
> |
q_i += *(idat.flucQ1); |
| 487 |
> |
} |
| 488 |
> |
|
| 489 |
|
if (idat.excluded) { |
| 490 |
|
*(idat.skippedCharge2) += q_i; |
| 491 |
|
} |
| 523 |
|
} |
| 524 |
|
|
| 525 |
|
if (j_is_Charge) { |
| 526 |
< |
q_j = data2.charge; |
| 526 |
> |
q_j = data2.fixedCharge; |
| 527 |
> |
|
| 528 |
> |
if (j_is_Fluctuating) |
| 529 |
> |
q_j += *(idat.flucQ2); |
| 530 |
> |
|
| 531 |
|
if (idat.excluded) { |
| 532 |
|
*(idat.skippedCharge1) += q_j; |
| 533 |
|
} |
| 565 |
|
duduz_j = V3Zero; |
| 566 |
|
} |
| 567 |
|
|
| 568 |
+ |
if (i_is_Fluctuating && j_is_Fluctuating) { |
| 569 |
+ |
J1 = Jij[idat.atypes]; |
| 570 |
+ |
J2 = Jij[make_pair(idat.atypes.second, idat.atypes.first)]; |
| 571 |
+ |
} |
| 572 |
+ |
|
| 573 |
|
epot = 0.0; |
| 574 |
|
dVdr = V3Zero; |
| 575 |
|
|
| 592 |
|
c2 = c1 * riji; |
| 593 |
|
} |
| 594 |
|
|
| 595 |
< |
preVal = *(idat.electroMult) * pre11_ * q_i * q_j; |
| 595 |
> |
preVal = *(idat.electroMult) * pre11_; |
| 596 |
|
|
| 597 |
|
if (summationMethod_ == esm_SHIFTED_POTENTIAL) { |
| 598 |
|
vterm = preVal * (c1 - c1c_); |
| 621 |
|
|
| 622 |
|
vterm = preVal * riji * erfcVal; |
| 623 |
|
dudr = - *(idat.sw) * preVal * c2; |
| 624 |
+ |
|
| 625 |
+ |
} |
| 626 |
+ |
|
| 627 |
+ |
vpair += vterm * q_i * q_j; |
| 628 |
+ |
epot += *(idat.sw) * vterm * q_i * q_j; |
| 629 |
+ |
dVdr += dudr * rhat * q_i * q_j; |
| 630 |
|
|
| 631 |
+ |
if (i_is_Fluctuating) { |
| 632 |
+ |
if (idat.excluded) { |
| 633 |
+ |
// vFluc1 is the difference between the direct coulomb integral |
| 634 |
+ |
// and the normal 1/r-like interaction between point charges. |
| 635 |
+ |
coulInt = J1->getValueAt( *(idat.rij) ); |
| 636 |
+ |
vFluc1 = coulInt - (*(idat.sw) * vterm); |
| 637 |
+ |
} else { |
| 638 |
+ |
vFluc1 = 0.0; |
| 639 |
+ |
} |
| 640 |
+ |
*(idat.dVdFQ1) += ( *(idat.sw) * vterm + vFluc1 ) * q_j; |
| 641 |
|
} |
| 642 |
|
|
| 643 |
< |
vpair += vterm; |
| 644 |
< |
epot += *(idat.sw) * vterm; |
| 645 |
< |
dVdr += dudr * rhat; |
| 643 |
> |
if (j_is_Fluctuating) { |
| 644 |
> |
if (idat.excluded) { |
| 645 |
> |
// vFluc2 is the difference between the direct coulomb integral |
| 646 |
> |
// and the normal 1/r-like interaction between point charges. |
| 647 |
> |
coulInt = J2->getValueAt( *(idat.rij) ); |
| 648 |
> |
vFluc2 = coulInt - (*(idat.sw) * vterm); |
| 649 |
> |
} else { |
| 650 |
> |
vFluc2 = 0.0; |
| 651 |
> |
} |
| 652 |
> |
*(idat.dVdFQ2) += ( *(idat.sw) * vterm + vFluc2 ) * q_i; |
| 653 |
> |
} |
| 654 |
> |
|
| 655 |
> |
|
| 656 |
|
} |
| 657 |
|
|
| 658 |
|
if (j_is_Dipole) { |
| 725 |
|
duduz_j += -preSw * pot_term * rhat; |
| 726 |
|
|
| 727 |
|
} |
| 728 |
+ |
if (i_is_Fluctuating) { |
| 729 |
+ |
*(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i; |
| 730 |
+ |
} |
| 731 |
|
} |
| 732 |
|
|
| 733 |
|
if (j_is_Quadrupole) { |
| 780 |
|
dudux_j += preSw * qxx_j * cx_j * rhatdot2; |
| 781 |
|
duduy_j += preSw * qyy_j * cy_j * rhatdot2; |
| 782 |
|
duduz_j += preSw * qzz_j * cz_j * rhatdot2; |
| 783 |
+ |
if (i_is_Fluctuating) { |
| 784 |
+ |
*(idat.dVdFQ1) += ( *(idat.sw) * vterm ) / q_i; |
| 785 |
+ |
} |
| 786 |
+ |
|
| 787 |
|
} |
| 788 |
|
} |
| 789 |
|
|
| 860 |
|
dVdr += preSw * (uz_i * c2ri - ct_i * rhat * sc2 * c3); |
| 861 |
|
duduz_i += preSw * pot_term * rhat; |
| 862 |
|
} |
| 863 |
+ |
|
| 864 |
+ |
if (j_is_Fluctuating) { |
| 865 |
+ |
*(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j; |
| 866 |
+ |
} |
| 867 |
+ |
|
| 868 |
|
} |
| 869 |
|
|
| 870 |
|
if (j_is_Dipole) { |
| 1015 |
|
dudux_i += preSw * qxx_i * cx_i * rhatdot2; |
| 1016 |
|
duduy_i += preSw * qyy_i * cy_i * rhatdot2; |
| 1017 |
|
duduz_i += preSw * qzz_i * cz_i * rhatdot2; |
| 1018 |
+ |
|
| 1019 |
+ |
if (j_is_Fluctuating) { |
| 1020 |
+ |
*(idat.dVdFQ2) += ( *(idat.sw) * vterm ) / q_j; |
| 1021 |
+ |
} |
| 1022 |
+ |
|
| 1023 |
|
} |
| 1024 |
|
} |
| 1025 |
|
|
| 1058 |
|
*(idat.t2) -= cross(uz_j, indirect_duduz_j); |
| 1059 |
|
} |
| 1060 |
|
|
| 937 |
– |
|
| 1061 |
|
return; |
| 1062 |
|
} |
| 1063 |
|
|
| 1064 |
|
void Electrostatic::calcSelfCorrection(SelfData &sdat) { |
| 1065 |
< |
RealType mu1, preVal, chg1, self; |
| 943 |
< |
|
| 1065 |
> |
RealType mu1, preVal, self; |
| 1066 |
|
if (!initialized_) initialize(); |
| 1067 |
|
|
| 1068 |
|
ElectrostaticAtomData data = ElectrostaticMap[sdat.atype]; |
| 1070 |
|
// logicals |
| 1071 |
|
bool i_is_Charge = data.is_Charge; |
| 1072 |
|
bool i_is_Dipole = data.is_Dipole; |
| 1073 |
+ |
bool i_is_Fluctuating = data.is_Fluctuating; |
| 1074 |
+ |
RealType chg1 = data.fixedCharge; |
| 1075 |
+ |
|
| 1076 |
+ |
if (i_is_Fluctuating) { |
| 1077 |
+ |
chg1 += *(sdat.flucQ); |
| 1078 |
+ |
// dVdFQ is really a force, so this is negative the derivative |
| 1079 |
+ |
*(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity; |
| 1080 |
+ |
} |
| 1081 |
|
|
| 1082 |
|
if (summationMethod_ == esm_REACTION_FIELD) { |
| 1083 |
|
if (i_is_Dipole) { |
| 1094 |
|
} |
| 1095 |
|
} else if (summationMethod_ == esm_SHIFTED_FORCE || summationMethod_ == esm_SHIFTED_POTENTIAL) { |
| 1096 |
|
if (i_is_Charge) { |
| 967 |
– |
chg1 = data.charge; |
| 1097 |
|
if (screeningMethod_ == DAMPED) { |
| 1098 |
|
self = - 0.5 * (c1c_ + alphaPi_) * chg1 * (chg1 + *(sdat.skippedCharge)) * pre11_; |
| 1099 |
|
} else { |
