Skip to content

Commit

Permalink
refactor: Clarify forces and derivatives (#1852)
Browse files Browse the repository at this point in the history
  • Loading branch information
@trisyoungs
trisyoungs committed May 9, 2024
1 parent 2584d89 commit 7d0de7e
Show file tree
Hide file tree
Showing 5 changed files with 109 additions and 103 deletions.
46 changes: 28 additions & 18 deletions src/classes/pairPotential.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ double PairPotential::analyticShortRangeEnergy(double r, PairPotential::ShortRan
// Apply the selected truncation scheme
if (truncation == PairPotential::ShiftedShortRangeTruncation)
{
energy += -(r - range_) * shortRangeForceAtCutoff_ - shortRangeEnergyAtCutoff_;
energy += (r - range_) * shortRangeForceAtCutoff_ - shortRangeEnergyAtCutoff_;
}

return energy;
Expand All @@ -212,11 +212,17 @@ double PairPotential::analyticShortRangeForce(double r, PairPotential::ShortRang
* Parameter 1 = Sigma
*/

// f = -48*epsilon*((sigma**12/x**13)-0.5*(sigma**6/x**7))
/*
* dU/dr = -48*epsilon*((sigma**12/x**13)-0.5*(sigma**6/x**7))
*
* signa**6 ( sigma**6 )
* = 48 * epsilon * -------- * ( - --------- + 0.5 ) * r**-1
* r**6 ( r**6 )
*/

auto sigmar = params[1] / r;
auto sigmar6 = pow(sigmar, 6.0);
force = 48.0 * params[0] * sigmar6 * (-sigmar6 + 0.5) / r;
force = -48.0 * params[0] * sigmar6 * (-sigmar6 + 0.5) / r;
}
break;
default:
Expand Down Expand Up @@ -369,13 +375,8 @@ double PairPotential::analyticEnergy(double r) const
if (r > range_)
return 0.0;

// Short-range potential
auto energy = analyticShortRangeEnergy(r);

// Coulomb contribution
energy += analyticCoulombEnergy(chargeI_ * chargeJ_, r);

return energy;
// Short-range potential and Coulomb contribution
return analyticShortRangeEnergy(r) + analyticCoulombEnergy(chargeI_ * chargeJ_, r);
}

// Return analytic potential at specified r, including Coulomb term from supplied charge product
Expand Down Expand Up @@ -405,7 +406,7 @@ double PairPotential::force(double r)
{
assert(r >= 0);

return dUFullInterpolation_.y(r, r * rDelta_);
return -dUFullInterpolation_.y(r, r * rDelta_);
}

// Return analytic force at specified r
Expand All @@ -414,11 +415,8 @@ double PairPotential::analyticForce(double r) const
if (r > range_)
return 0.0;

// Short-range potential
double force = analyticShortRangeForce(r);

// Coulomb contribution
force += analyticCoulombForce(chargeI_ * chargeJ_, r);
// Short-range potential and Coulomb contribution
auto force = analyticShortRangeForce(r) + analyticCoulombForce(chargeI_ * chargeJ_, r);

return force;
}
Expand All @@ -436,11 +434,23 @@ double PairPotential::analyticForce(double qiqj, double r, double elecScale, dou
// Return analytic coulomb force of specified charges
double PairPotential::analyticCoulombForce(double qiqj, double r, PairPotential::CoulombTruncationScheme truncation) const
{
/*
* Derivative of Coulomb's Law is:
*
* dU q(i)*q(j) q(i)*q(j)
* dU/dr = -- --------- = - ---------
* dr r r*r
*
* q(i)*q(j)
* The force is -(dU/dr) = ---------
* r*r
*/

// Calculate based on truncation scheme
if (truncation == PairPotential::NoCoulombTruncation)
return -COULCONVERT * qiqj / (r * r);
return COULCONVERT * qiqj / (r * r);
else if (truncation == PairPotential::ShiftedCoulombTruncation)
return -COULCONVERT * qiqj * (1.0 / (r * r) - 1.0 / (range_ * range_));
return COULCONVERT * qiqj * (1.0 / (r * r) - 1.0 / (range_ * range_));

return 0.0;
}
Expand Down
4 changes: 2 additions & 2 deletions src/io/export/pairPotential.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,10 @@ bool PairPotentialExportFileFormat::exportBlock(LineParser &parser, PairPotentia

// Write header comment
if (!parser.writeLineF("#{:9} {:12} {:12} {:12} {:12} {:12} {:12}\n", "", "Full", "Derivative", "Original",
"Additional", "Exact(Orig)", "Exact(Deriv)"))
"Additional", "Exact(Orig)", "Exact(Force)"))
return false;
if (!parser.writeLineF("#{:9} {:12} {:12} {:12} {:12} {:12} {:12}\n", "r(Angs)", "U(kJ/mol)", "dU(kJ/mol/Ang)",
"U(kJ/mol)", "U(kJ/mol)", "U(kJ/mol)", "dU(kJ/mol/Ang)"))
"U(kJ/mol)", "U(kJ/mol)", "U(kJ/mol)", "F(kJ/mol/Ang)"))
return false;

for (auto n = 0; n < nPoints; ++n)
Expand Down
48 changes: 24 additions & 24 deletions src/kernels/force.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,59 +30,59 @@ ForceKernel::ForceKernel(const Box *box, const ProcessPool &procPool, const Pote
// Calculate PairPotential forces between Atoms provided
void ForceKernel::forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f) const
{
auto force = j.r() - i.r();
auto distanceSq = force.magnitudeSq();
auto vij = j.r() - i.r();
auto distanceSq = vij.magnitudeSq();
if (distanceSq > cutoffDistanceSquared_)
return;
auto r = sqrt(distanceSq);
force /= r;
force *= potentialMap_.force(i, j, r);
f[indexI] += force;
f[indexJ] -= force;
vij /= r;
vij *= potentialMap_.force(i, j, r);
f[indexI] -= vij;
f[indexJ] += vij;
}

// Calculate inter-particle forces between Atoms provided, scaling electrostatic and van der Waals components
void ForceKernel::forcesWithoutMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f, double elecScale,
double vdwScale) const
{
auto force = j.r() - i.r();
auto distanceSq = force.magnitudeSq();
auto vij = j.r() - i.r();
auto distanceSq = vij.magnitudeSq();
if (distanceSq > cutoffDistanceSquared_)
return;
auto r = sqrt(distanceSq);
force /= r;
force *= potentialMap_.force(i, j, r, elecScale, vdwScale);
f[indexI] += force;
f[indexJ] -= force;
vij /= r;
vij *= potentialMap_.force(i, j, r, elecScale, vdwScale);
f[indexI] -= vij;
f[indexJ] += vij;
}

// Calculate PairPotential forces between Atoms provided
void ForceKernel::forcesWithMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f) const
{
auto force = box_->minimumVector(i.r(), j.r());
auto distanceSq = force.magnitudeSq();
auto vij = box_->minimumVector(i.r(), j.r());
auto distanceSq = vij.magnitudeSq();
if (distanceSq > cutoffDistanceSquared_)
return;
auto r = sqrt(distanceSq);
force /= r;
force *= potentialMap_.force(i, j, r);
f[indexI] += force;
f[indexJ] -= force;
vij /= r;
vij *= potentialMap_.force(i, j, r);
f[indexI] -= vij;
f[indexJ] += vij;
}

// Calculate inter-particle forces between Atoms provided, scaling electrostatic and van der Waals components
void ForceKernel::forcesWithMim(const Atom &i, int indexI, const Atom &j, int indexJ, ForceVector &f, double elecScale,
double vdwScale) const
{
auto force = box_->minimumVector(i.r(), j.r());
auto distanceSq = force.magnitudeSq();
auto vij = box_->minimumVector(i.r(), j.r());
auto distanceSq = vij.magnitudeSq();
if (distanceSq > cutoffDistanceSquared_)
return;
auto r = sqrt(distanceSq);
force /= r;
force *= potentialMap_.force(i, j, r, elecScale, vdwScale);
f[indexI] += force;
f[indexJ] -= force;
vij /= r;
vij *= potentialMap_.force(i, j, r, elecScale, vdwScale);
f[indexI] -= vij;
f[indexJ] += vij;
}

/*
Expand Down
4 changes: 2 additions & 2 deletions src/modules/forces/functions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,8 +115,8 @@ void ForcesModule::totalForces(const ProcessPool &procPool, const Species *sp, c
vecij *= potentialMap.force(&i, &j, r, elec14, vdw14);

auto &fLocal = combinableUnbound.local();
fLocal[indexI] += vecij;
fLocal[indexJ] -= vecij;
fLocal[indexI] -= vecij;
fLocal[indexJ] += vecij;
};

// Calculate pairwise forces between atoms
Expand Down

1 comment on commit 7d0de7e

@github-actions
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Performance Alert ⚠️

Possible performance regression was detected for benchmark.
Benchmark result of this commit is worse than the previous benchmark result exceeding threshold 2.

Benchmark suite Current: 7d0de7e Previous: 2584d89 Ratio
BM_HistogramBinning_2d/16777216 51.13413653322156 ns/iter 20.195253776281803 ns/iter 2.53

This comment was automatically generated by workflow using github-action-benchmark.

CC: @disorderedmaterials/dissolve-devs

Please sign in to comment.