321 |
|
electrostaticAtomData.hardness = fqa.getHardness(); |
322 |
|
electrostaticAtomData.slaterN = fqa.getSlaterN(); |
323 |
|
electrostaticAtomData.slaterZeta = fqa.getSlaterZeta(); |
324 |
+ |
} else { |
325 |
+ |
electrostaticAtomData.is_Fluctuating = false; |
326 |
|
} |
327 |
|
|
328 |
|
pair<map<int,AtomType*>::iterator,bool> ret; |
366 |
|
rval = RealType(i) * dr; |
367 |
|
rvals.push_back(rval); |
368 |
|
J1vals.push_back( sSTOCoulInt( a, b, m, n, rval * PhysicalConstants::angstromsToBohr ) ); |
369 |
+ |
// may not be necessary if Slater coulomb integral is symmetric |
370 |
|
J2vals.push_back( sSTOCoulInt( b, a, n, m, rval * PhysicalConstants::angstromsToBohr ) ); |
371 |
|
} |
372 |
|
|
450 |
|
Vector3d indirect_duduz_i(V3Zero), indirect_duduz_j(V3Zero); |
451 |
|
|
452 |
|
pair<RealType, RealType> res; |
453 |
+ |
|
454 |
+ |
// splines for coulomb integrals |
455 |
+ |
CubicSpline* J1; |
456 |
+ |
CubicSpline* J2; |
457 |
|
|
458 |
|
if (!initialized_) initialize(); |
459 |
|
|
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.fixedCharge; |
484 |
+ |
|
485 |
+ |
if (i_is_Fluctuating) { |
486 |
+ |
q_i += *(idat.flucQ1); |
487 |
+ |
} |
488 |
+ |
|
489 |
|
if (idat.excluded) { |
490 |
|
*(idat.skippedCharge2) += q_i; |
491 |
|
} |
524 |
|
|
525 |
|
if (j_is_Charge) { |
526 |
|
q_j = data2.fixedCharge; |
527 |
+ |
|
528 |
+ |
if (i_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 |
|
|
621 |
|
|
622 |
|
vterm = preVal * riji * erfcVal; |
623 |
|
dudr = - *(idat.sw) * preVal * c2; |
624 |
+ |
|
625 |
+ |
} |
626 |
|
|
627 |
+ |
|
628 |
+ |
if (i_is_Fluctuating) { |
629 |
+ |
if (!idat.excluded) |
630 |
+ |
*(idat.dVdFQ1) += *(idat.sw) * vterm / q_i; |
631 |
+ |
else { |
632 |
+ |
res = J1->getValueAndDerivativeAt( *(idat.rij) ); |
633 |
+ |
*(idat.dVdFQ1) += pre11_ * res.first * q_j; |
634 |
+ |
} |
635 |
|
} |
636 |
+ |
if (j_is_Fluctuating) { |
637 |
+ |
if (!idat.excluded) |
638 |
+ |
*(idat.dVdFQ2) += *(idat.sw) * vterm / q_j; |
639 |
+ |
else { |
640 |
+ |
res = J2->getValueAndDerivativeAt( *(idat.rij) ); |
641 |
+ |
*(idat.dVdFQ2) += pre11_ * res.first * q_i; |
642 |
+ |
} |
643 |
+ |
} |
644 |
|
|
645 |
|
vpair += vterm; |
646 |
|
epot += *(idat.sw) * vterm; |