64 pairPotentialProperties_(pairPotentialProperties),
65 rCut_(
readScalar(pairPotentialProperties_.lookup(
"rCut"))),
66 rCutSqr_(rCut_*rCut_),
67 rMin_(
readScalar(pairPotentialProperties_.lookup(
"rMin"))),
68 dr_(
readScalar(pairPotentialProperties_.lookup(
"dr"))),
72 writeTables_(
Switch(pairPotentialProperties_.lookup(
"writeTables")))
80 label N = label((rCut_ - rMin_)/dr_) + 1;
82 forceLookup_.setSize(N);
84 energyLookup_.setSize(N);
88 energyLookup_[
k] = scaledEnergy(
k*dr_ + rMin_);
90 forceLookup_[
k] = -energyDerivative((
k*dr_ + rMin_),
true);
97 scalar k_rIJ = (r - rMin_)/dr_;
99 label
k = label(k_rIJ);
104 <<
"r less than rMin in pair potential " << name_ <<
nl
109 (k_rIJ -
k)*forceLookup_[k+1]
110 + (k + 1 - k_rIJ)*forceLookup_[
k];
123 forceTab[
k].first() = rMin_ +
k*dr_;
125 forceTab[
k].second() = forceLookup_[
k];
134 scalar k_rIJ = (r - rMin_)/dr_;
136 label
k = label(k_rIJ);
141 <<
"r less than rMin in pair potential " << name_ <<
nl
146 (k_rIJ -
k)*energyLookup_[k+1]
147 + (k + 1 - k_rIJ)*energyLookup_[
k];
160 energyTab[
k].first() = rMin_ +
k*dr_;
162 energyTab[
k].second() = energyLookup_[
k];
171 scalar
e = unscaledEnergy(r);
182 const bool scaledEnergyDerivative
194 if (scaledEnergyDerivative)
196 Ea = scaledEnergy(ra);
197 Ef = scaledEnergy(rf);
198 Eb = scaledEnergy(rb);
202 Ea = unscaledEnergy(ra);
203 Ef = unscaledEnergy(rf);
204 Eb = unscaledEnergy(rb);
207 scalar denominator = (ra - rf)*(ra - rb)*(rf - rb);
211 rb*rb*(Ea - Ef) + ra*ra*(Ef - Eb) + rf*rf*(Eb - Ea)
216 rb*(Ef - Ea) + rf*(Ea - Eb) + ra*(Eb - Ef)
219 return a1 + 2.0*a2*r;
225 pairPotentialProperties_ = pairPotentialProperties;