| 743 |
|
} |
| 744 |
|
|
| 745 |
|
void Electrostatic::calcForce(InteractionData &idat) { |
| 746 |
– |
|
| 747 |
– |
RealType C_a, C_b; // Charges |
| 748 |
– |
Vector3d D_a, D_b; // Dipoles (space-fixed) |
| 749 |
– |
Mat3x3d Q_a, Q_b; // Quadrupoles (space-fixed) |
| 746 |
|
|
| 747 |
< |
RealType ri; // Distance utility scalar |
| 748 |
< |
RealType rdDa, rdDb; // Dipole utility scalars |
| 749 |
< |
Vector3d rxDa, rxDb; // Dipole utility vectors |
| 750 |
< |
RealType rdQar, rdQbr, trQa, trQb; // Quadrupole utility scalars |
| 755 |
< |
Vector3d Qar, Qbr, rQa, rQb, rxQar, rxQbr; // Quadrupole utility vectors |
| 756 |
< |
RealType pref; |
| 757 |
< |
|
| 758 |
< |
RealType DadDb, trQaQb, DadQbr, DbdQar; // Cross-interaction scalars |
| 759 |
< |
RealType rQaQbr; |
| 760 |
< |
Vector3d DaxDb, DadQb, DbdQa, DaxQbr, DbxQar; // Cross-interaction vectors |
| 761 |
< |
Vector3d rQaQb, QaQbr, QaxQb, rQaxQbr; |
| 762 |
< |
Mat3x3d QaQb; // Cross-interaction matrices |
| 747 |
> |
if (!initialized_) initialize(); |
| 748 |
> |
|
| 749 |
> |
data1 = ElectrostaticMap[idat.atypes.first]; |
| 750 |
> |
data2 = ElectrostaticMap[idat.atypes.second]; |
| 751 |
|
|
| 752 |
< |
RealType U(0.0); // Potential |
| 753 |
< |
Vector3d F(0.0); // Force |
| 754 |
< |
Vector3d Ta(0.0); // Torque on site a |
| 755 |
< |
Vector3d Tb(0.0); // Torque on site b |
| 756 |
< |
Vector3d Ea(0.0); // Electric field at site a |
| 757 |
< |
Vector3d Eb(0.0); // Electric field at site b |
| 758 |
< |
RealType dUdCa(0.0); // fluctuating charge force at site a |
| 759 |
< |
RealType dUdCb(0.0); // fluctuating charge force at site a |
| 752 |
> |
U = 0.0; // Potential |
| 753 |
> |
F.zero(); // Force |
| 754 |
> |
Ta.zero(); // Torque on site a |
| 755 |
> |
Tb.zero(); // Torque on site b |
| 756 |
> |
Ea.zero(); // Electric field at site a |
| 757 |
> |
Eb.zero(); // Electric field at site b |
| 758 |
> |
dUdCa = 0.0; // fluctuating charge force at site a |
| 759 |
> |
dUdCb = 0.0; // fluctuating charge force at site a |
| 760 |
|
|
| 761 |
|
// Indirect interactions mediated by the reaction field. |
| 762 |
< |
RealType indirect_Pot(0.0); // Potential |
| 763 |
< |
Vector3d indirect_F(0.0); // Force |
| 764 |
< |
Vector3d indirect_Ta(0.0); // Torque on site a |
| 765 |
< |
Vector3d indirect_Tb(0.0); // Torque on site b |
| 762 |
> |
indirect_Pot = 0.0; // Potential |
| 763 |
> |
indirect_F.zero(); // Force |
| 764 |
> |
indirect_Ta.zero(); // Torque on site a |
| 765 |
> |
indirect_Tb.zero(); // Torque on site b |
| 766 |
|
|
| 767 |
|
// Excluded potential that is still computed for fluctuating charges |
| 768 |
< |
RealType excluded_Pot(0.0); |
| 768 |
> |
excluded_Pot= 0.0; |
| 769 |
|
|
| 782 |
– |
RealType rfContrib, coulInt; |
| 783 |
– |
|
| 784 |
– |
// spline for coulomb integral |
| 785 |
– |
CubicSpline* J; |
| 770 |
|
|
| 787 |
– |
if (!initialized_) initialize(); |
| 788 |
– |
|
| 789 |
– |
ElectrostaticAtomData data1 = ElectrostaticMap[idat.atypes.first]; |
| 790 |
– |
ElectrostaticAtomData data2 = ElectrostaticMap[idat.atypes.second]; |
| 791 |
– |
|
| 771 |
|
// some variables we'll need independent of electrostatic type: |
| 772 |
|
|
| 773 |
|
ri = 1.0 / *(idat.rij); |
| 774 |
< |
Vector3d rhat = *(idat.d) * ri; |
| 774 |
> |
rhat = *(idat.d) * ri; |
| 775 |
|
|
| 776 |
|
// logicals |
| 777 |
|
|
| 778 |
< |
bool a_is_Charge = data1.is_Charge; |
| 779 |
< |
bool a_is_Dipole = data1.is_Dipole; |
| 780 |
< |
bool a_is_Quadrupole = data1.is_Quadrupole; |
| 781 |
< |
bool a_is_Fluctuating = data1.is_Fluctuating; |
| 778 |
> |
a_is_Charge = data1.is_Charge; |
| 779 |
> |
a_is_Dipole = data1.is_Dipole; |
| 780 |
> |
a_is_Quadrupole = data1.is_Quadrupole; |
| 781 |
> |
a_is_Fluctuating = data1.is_Fluctuating; |
| 782 |
|
|
| 783 |
< |
bool b_is_Charge = data2.is_Charge; |
| 784 |
< |
bool b_is_Dipole = data2.is_Dipole; |
| 785 |
< |
bool b_is_Quadrupole = data2.is_Quadrupole; |
| 786 |
< |
bool b_is_Fluctuating = data2.is_Fluctuating; |
| 783 |
> |
b_is_Charge = data2.is_Charge; |
| 784 |
> |
b_is_Dipole = data2.is_Dipole; |
| 785 |
> |
b_is_Quadrupole = data2.is_Quadrupole; |
| 786 |
> |
b_is_Fluctuating = data2.is_Fluctuating; |
| 787 |
|
|
| 788 |
|
// Obtain all of the required radial function values from the |
| 789 |
|
// spline structures: |
