| 1555 |
|
mMin = 1; |
| 1556 |
|
} |
| 1557 |
|
pot += kPot; |
| 1558 |
+ |
} |
| 1559 |
+ |
|
| 1560 |
+ |
void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded, |
| 1561 |
+ |
RealType &spot1, RealType &spot2) { |
| 1562 |
+ |
|
| 1563 |
+ |
if (!initialized_) { |
| 1564 |
+ |
cerr << "initializing\n"; |
| 1565 |
+ |
initialize(); |
| 1566 |
+ |
cerr << "done\n"; |
| 1567 |
+ |
} |
| 1568 |
+ |
|
| 1569 |
+ |
const RealType mPoleConverter = 0.20819434; |
| 1570 |
+ |
|
| 1571 |
+ |
AtomType* atype1 = a1->getAtomType(); |
| 1572 |
+ |
AtomType* atype2 = a2->getAtomType(); |
| 1573 |
+ |
int atid1 = atype1->getIdent(); |
| 1574 |
+ |
int atid2 = atype2->getIdent(); |
| 1575 |
+ |
data1 = ElectrostaticMap[Etids[atid1]]; |
| 1576 |
+ |
data2 = ElectrostaticMap[Etids[atid2]]; |
| 1577 |
+ |
|
| 1578 |
+ |
Pa = 0.0; // Site potential at site a |
| 1579 |
+ |
Pb = 0.0; // Site potential at site b |
| 1580 |
+ |
|
| 1581 |
+ |
Vector3d d = a2->getPos() - a1->getPos(); |
| 1582 |
+ |
info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d); |
| 1583 |
+ |
RealType rij = d.length(); |
| 1584 |
+ |
// some variables we'll need independent of electrostatic type: |
| 1585 |
+ |
|
| 1586 |
+ |
RealType ri = 1.0 / rij; |
| 1587 |
+ |
rhat = d * ri; |
| 1588 |
+ |
|
| 1589 |
+ |
|
| 1590 |
+ |
if ((rij >= cutoffRadius_) || excluded) { |
| 1591 |
+ |
spot1 = 0.0; |
| 1592 |
+ |
spot2 = 0.0; |
| 1593 |
+ |
return; |
| 1594 |
+ |
} |
| 1595 |
+ |
|
| 1596 |
+ |
// logicals |
| 1597 |
+ |
|
| 1598 |
+ |
a_is_Charge = data1.is_Charge; |
| 1599 |
+ |
a_is_Dipole = data1.is_Dipole; |
| 1600 |
+ |
a_is_Quadrupole = data1.is_Quadrupole; |
| 1601 |
+ |
a_is_Fluctuating = data1.is_Fluctuating; |
| 1602 |
+ |
|
| 1603 |
+ |
b_is_Charge = data2.is_Charge; |
| 1604 |
+ |
b_is_Dipole = data2.is_Dipole; |
| 1605 |
+ |
b_is_Quadrupole = data2.is_Quadrupole; |
| 1606 |
+ |
b_is_Fluctuating = data2.is_Fluctuating; |
| 1607 |
+ |
|
| 1608 |
+ |
// Obtain all of the required radial function values from the |
| 1609 |
+ |
// spline structures: |
| 1610 |
+ |
|
| 1611 |
+ |
|
| 1612 |
+ |
if (a_is_Charge || b_is_Charge) { |
| 1613 |
+ |
v01 = v01s->getValueAt(rij); |
| 1614 |
+ |
} |
| 1615 |
+ |
if (a_is_Dipole || b_is_Dipole) { |
| 1616 |
+ |
v11 = v11s->getValueAt(rij); |
| 1617 |
+ |
v11or = ri * v11; |
| 1618 |
+ |
} |
| 1619 |
+ |
if (a_is_Quadrupole || b_is_Quadrupole) { |
| 1620 |
+ |
v21 = v21s->getValueAt(rij); |
| 1621 |
+ |
v22 = v22s->getValueAt(rij); |
| 1622 |
+ |
v22or = ri * v22; |
| 1623 |
+ |
} |
| 1624 |
+ |
|
| 1625 |
+ |
if (a_is_Charge) { |
| 1626 |
+ |
C_a = data1.fixedCharge; |
| 1627 |
+ |
|
| 1628 |
+ |
if (a_is_Fluctuating) { |
| 1629 |
+ |
C_a += a1->getFlucQPos(); |
| 1630 |
+ |
} |
| 1631 |
+ |
|
| 1632 |
+ |
Pb += C_a * pre11_ * v01; |
| 1633 |
+ |
} |
| 1634 |
+ |
|
| 1635 |
+ |
if (a_is_Dipole) { |
| 1636 |
+ |
D_a = a1->getDipole() * mPoleConverter; |
| 1637 |
+ |
rdDa = dot(rhat, D_a); |
| 1638 |
+ |
Pb += pre12_ * v11 * rdDa; |
| 1639 |
+ |
} |
| 1640 |
+ |
|
| 1641 |
+ |
if (a_is_Quadrupole) { |
| 1642 |
+ |
Q_a = a1->getQuadrupole() * mPoleConverter; |
| 1643 |
+ |
trQa = Q_a.trace(); |
| 1644 |
+ |
Qar = Q_a * rhat; |
| 1645 |
+ |
rdQar = dot(rhat, Qar); |
| 1646 |
+ |
Pb += pre14_ * (v21 * trQa + v22 * rdQar); |
| 1647 |
+ |
} |
| 1648 |
+ |
|
| 1649 |
+ |
if (b_is_Charge) { |
| 1650 |
+ |
C_b = data2.fixedCharge; |
| 1651 |
+ |
|
| 1652 |
+ |
if (b_is_Fluctuating) |
| 1653 |
+ |
C_b += a2->getFlucQPos(); |
| 1654 |
+ |
|
| 1655 |
+ |
Pa += C_b * pre11_ * v01; |
| 1656 |
+ |
} |
| 1657 |
+ |
|
| 1658 |
+ |
if (b_is_Dipole) { |
| 1659 |
+ |
D_b = a2->getDipole() * mPoleConverter; |
| 1660 |
+ |
rdDb = dot(rhat, D_b); |
| 1661 |
+ |
Pa += pre12_ * v11 * rdDb; |
| 1662 |
+ |
} |
| 1663 |
+ |
|
| 1664 |
+ |
if (b_is_Quadrupole) { |
| 1665 |
+ |
Q_a = a2->getQuadrupole() * mPoleConverter; |
| 1666 |
+ |
trQb = Q_b.trace(); |
| 1667 |
+ |
Qbr = Q_b * rhat; |
| 1668 |
+ |
rdQbr = dot(rhat, Qbr); |
| 1669 |
+ |
Pa += pre14_ * (v21 * trQb + v22 * rdQbr); |
| 1670 |
+ |
} |
| 1671 |
+ |
|
| 1672 |
+ |
spot1 = Pa; |
| 1673 |
+ |
spot2 = Pb; |
| 1674 |
|
} |
| 1675 |
|
} |
