Skip to content

Commit 9451266

Browse files
jgreener64jgreener64
authored
Merge pull request #232 from Immi000/lambda_factor
Add missing λ-factor in force and potential for LennardJonesSoftCoreBeutler
2 parents 377a405 + 3ea323b commit 9451266

File tree

3 files changed

+50
-50
lines changed

3 files changed

+50
-50
lines changed

‎src/interactions/coulomb.jl‎

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -118,15 +118,15 @@ the appearing and disappearing of atoms.
118118
See [Beutler et al. 1994](https://doi.org/10.1016/0009-2614(94)00397-1).
119119
The potential energy is defined as
120120
```math
121-
V(r_{ij}) = \frac{1}{4\pi\epsilon_0} \frac{q_iq_j}{r_Q^{1/6}}
121+
V(r_{ij}) = \lambda\frac{1}{4\pi\epsilon_0} \frac{q_iq_j}{r_Q^{1/6}}
122122
```
123123
and the force on each atom by
124124
```math
125-
\vec{F}_i = \frac{1}{4\pi\epsilon_0} \frac{q_iq_jr_{ij}^5}{r_Q^{7/6}}\frac{\vec{r_{ij}}}{r_{ij}}
125+
\vec{F}_i = \lambda\frac{1}{4\pi\epsilon_0} \frac{q_iq_jr_{ij}^5}{r_Q^{7/6}}\frac{\vec{r_{ij}}}{r_{ij}}
126126
```
127127
where
128128
```math
129-
r_{Q} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}})+r_{ij}^6)
129+
r_{Q} = \left(\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}\right)+r_{ij}^6
130130
```
131131
and
132132
```math
@@ -194,7 +194,7 @@ end
194194
σ6 = inter.σ_mixing(atom_i, atom_j)^6
195195
ϵ = inter.ϵ_mixing(atom_i, atom_j)
196196
C6 = 4 * ϵ * σ6
197-
params = (ke, qi, qj, C6 * σ6, C6, inter.σ6_fac)
197+
params = (ke, qi, qj, C6 * σ6, C6, inter.σ6_fac, inter.λ)
198198

199199
f = force_cutoff(cutoff, inter, r, params)
200200
fdr = (f / r) * dr
@@ -205,10 +205,10 @@ end
205205
end
206206
end
207207

208-
function pairwise_force(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac))
208+
function pairwise_force(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac, λ))
209209
r3 = r^3
210210
R = ((σ6_fac*(C12/C6))+(r3*r3))*sqrt(cbrt(((σ6_fac*(C12/C6))+(r3*r3))))
211-
return ke * ((qi*qj)/R) * (r3*r*r)
211+
return λ * ke * ((qi*qj)/R) * (r3*r*r)
212212
end
213213

214214
@inline function potential_energy(inter::CoulombSoftCoreBeutler,
@@ -225,7 +225,7 @@ end
225225
σ6 = inter.σ_mixing(atom_i, atom_j)^6
226226
ϵ = inter.ϵ_mixing(atom_i, atom_j)
227227
C6 = 4 * ϵ * σ6
228-
params = (ke, qi, qj, C6 *σ6, C6, inter.σ6_fac)
228+
params = (ke, qi, qj, C6 *σ6, C6, inter.σ6_fac, inter.λ)
229229

230230
pe = pe_cutoff(cutoff, inter, r, params)
231231
if special
@@ -235,9 +235,9 @@ end
235235
end
236236
end
237237

238-
function pairwise_pe(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac))
238+
function pairwise_pe(::CoulombSoftCoreBeutler, r, (ke, qi, qj, C12, C6, σ6_fac, λ))
239239
R = sqrt(cbrt((σ6_fac*(C12/C6))+r^6))
240-
return ke * ((qi * qj)/R)
240+
return λ * ke * ((qi * qj)/R)
241241
end
242242

243243
@doc raw"""
@@ -321,7 +321,7 @@ end
321321
cutoff = inter.cutoff
322322
ke = inter.coulomb_const
323323
qi, qj = atom_i.charge, atom_j.charge
324-
params = (ke, qi, qj, inter.σQ, inter.σ6_fac)
324+
params = (ke, qi, qj, inter.σQ, inter.σ6_fac, inter.λ)
325325

326326
f = force_cutoff(cutoff, inter, r, params)
327327
fdr = (f / r) * dr
@@ -332,13 +332,13 @@ end
332332
end
333333
end
334334

335-
function pairwise_force(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac))
335+
function pairwise_force(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac, λ))
336336
qij = qi * qj
337337
R = σ6_fac*(oneunit(r)+(σQ*abs(qij)))
338338
if r >= R
339-
return ke * (qij/(r^2))
339+
return λ * ke * (qij/(r^2))
340340
elseif r < R
341-
return ke * (-(((2*qij)/(R^3)) * r) + ((3*qij)/(R^2)))
341+
return λ * ke * (-(((2*qij)/(R^3)) * r) + ((3*qij)/(R^2)))
342342
end
343343
end
344344

@@ -353,7 +353,7 @@ end
353353
cutoff = inter.cutoff
354354
ke = inter.coulomb_const
355355
qi, qj = atom_i.charge, atom_j.charge
356-
params = (ke, qi, qj, inter.σQ, inter.σ6_fac)
356+
params = (ke, qi, qj, inter.σQ, inter.σ6_fac, inter.λ)
357357

358358
pe = pe_cutoff(cutoff, inter, r, params)
359359
if special
@@ -363,13 +363,13 @@ end
363363
end
364364
end
365365

366-
function pairwise_pe(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac))
366+
function pairwise_pe(::CoulombSoftCoreGapsys, r, (ke, qi, qj, σQ, σ6_fac, λ))
367367
qij = qi * qj
368368
R = σ6_fac*(oneunit(r)+(σQ*abs(qij)))
369369
if r >= R
370-
return ke * (qij/r)
370+
return λ * ke * (qij/r)
371371
elseif r < R
372-
return ke * (((qij/(R^3))*(r^2))-(((3*qij)/(R^2))*r)+((3*qij)/R))
372+
return λ * ke * (((qij/(R^3))*(r^2))-(((3*qij)/(R^2))*r)+((3*qij)/R))
373373
end
374374
end
375375

‎src/interactions/lennard_jones.jl‎

Lines changed: 17 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -149,15 +149,15 @@ the appearing and disappearing of atoms.
149149
See [Beutler et al. 1994](https://doi.org/10.1016/0009-2614(94)00397-1).
150150
The potential energy is defined as
151151
```math
152-
V(r_{ij}) = \frac{C^{(12)}}{r_{LJ}^{12}} - \frac{C^{(6)}}{r_{LJ}^{6}}
152+
V(r_{ij}) = \lambda\left(\frac{C^{(12)}}{r_{LJ}^{12}} - \frac{C^{(6)}}{r_{LJ}^{6}}\right
153153
```
154154
and the force on each atom by
155155
```math
156-
\vec{F}_i = ((\frac{12C^{(12)}}{r_{LJ}^{13}} - \frac{6C^{(6)}}{r_{LJ}^7})(\frac{r_{ij}}{r_{LJ}})^5) \frac{\vec{r_{ij}}}{r_{ij}}
156+
\vec{F}_i = \lambda\left(\left(\frac{12C^{(12)}}{r_{LJ}^{13}} - \frac{6C^{(6)}}{r_{LJ}^7}\right)\left(\frac{r_{ij}}{r_{LJ}}\right)^5\right) \frac{\vec{r_{ij}}}{r_{ij}}
157157
```
158158
where
159159
```math
160-
r_{LJ} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r^6)^{1/6}
160+
r_{LJ} = \left(\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r^6\right)^{1/6}
161161
```
162162
and
163163
```math
@@ -228,7 +228,7 @@ end
228228
cutoff = inter.cutoff
229229
r = norm(dr)
230230
C6 = 4 * ϵ * σ6
231-
params = (C6 * σ6, C6, inter.σ6_fac)
231+
params = (C6 * σ6, C6, inter.σ6_fac, inter.λ)
232232

233233
f = force_cutoff(cutoff, inter, r, params)
234234
fdr = (f / r) * dr
@@ -239,10 +239,10 @@ end
239239
end
240240
end
241241

242-
function pairwise_force(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac))
242+
function pairwise_force(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac, λ))
243243
R = sqrt(cbrt((σ6_fac*(C12/C6))+r^6))
244244
R6 = R^6
245-
return (((12*C12)/(R6*R6*R)) - ((6*C6)/(R6*R)))*((r/R)^5)
245+
return λ*(((12*C12)/(R6*R6*R)) - ((6*C6)/(R6*R)))*((r/R)^5)
246246
end
247247

248248
@inline function potential_energy(inter::LennardJonesSoftCoreBeutler,
@@ -261,7 +261,7 @@ end
261261
cutoff = inter.cutoff
262262
r = norm(dr)
263263
C6 = 4 * ϵ * σ6
264-
params = (C6 * σ6, C6, inter.σ6_fac)
264+
params = (C6 * σ6, C6, inter.σ6_fac, inter.λ)
265265

266266
pe = pe_cutoff(cutoff, inter, r, params)
267267
if special
@@ -271,9 +271,9 @@ end
271271
end
272272
end
273273

274-
function pairwise_pe(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac))
274+
function pairwise_pe(::LennardJonesSoftCoreBeutler, r, (C12, C6, σ6_fac, λ))
275275
R6 = (σ6_fac*(C12/C6))+r^6
276-
return ((C12/(R6*R6)) - (C6/(R6)))
276+
return λ*((C12/(R6*R6)) - (C6/(R6)))
277277
end
278278

279279
@doc raw"""
@@ -368,7 +368,7 @@ end
368368
cutoff = inter.cutoff
369369
r = norm(dr)
370370
C6 = 4 * ϵ * σ6
371-
params = (C6 * σ6, C6)
371+
params = (C6 * σ6, C6, inter.λ)
372372

373373
f = force_cutoff(cutoff, inter, r, params)
374374
fdr = (f / r) * dr
@@ -379,16 +379,16 @@ end
379379
end
380380
end
381381

382-
function pairwise_force(inter::LennardJonesSoftCoreGapsys, r, (C12, C6))
382+
function pairwise_force(inter::LennardJonesSoftCoreGapsys, r, (C12, C6, λ))
383383
R = inter.α*sqrt(cbrt((26*(C12/C6)*(1-inter.λ)/7)))
384384
r6 = r^6
385385
invR = inv(R)
386386
invR2 = invR^2
387387
invR6 = invR^6
388388
if r >= R
389-
return (((12*C12)/(r6*r6*r))-((6*C6)/(r6*r)))
389+
return λ * (((12*C12)/(r6*r6*r))-((6*C6)/(r6*r)))
390390
elseif r < R
391-
return (((-156*C12*(invR6*invR6*invR2)) + (42*C6*(invR2*invR6)))*r +
391+
return λ * (((-156*C12*(invR6*invR6*invR2)) + (42*C6*(invR2*invR6)))*r +
392392
(168*C12*(invR6*invR6*invR)) - (48*C6*(invR6*invR)))
393393
end
394394
end
@@ -409,7 +409,7 @@ end
409409
cutoff = inter.cutoff
410410
r = norm(dr)
411411
C6 = 4 * ϵ * σ6
412-
params = (C6 * σ6, C6)
412+
params = (C6 * σ6, C6, inter.λ)
413413

414414
pe = pe_cutoff(cutoff, inter, r, params)
415415
if special
@@ -419,16 +419,16 @@ end
419419
end
420420
end
421421

422-
function pairwise_pe(inter::LennardJonesSoftCoreGapsys, r, (C12, C6))
422+
function pairwise_pe(inter::LennardJonesSoftCoreGapsys, r, (C12, C6, λ))
423423
R = inter.α*sqrt(cbrt((26*(C12/C6)*(1-inter.λ)/7)))
424424
r6 = r^6
425425
invR = inv(R)
426426
invR2 = invR^2
427427
invR6 = invR^6
428428
if r >= R
429-
return (C12/(r6*r6))-(C6/(r6))
429+
return λ * (C12/(r6*r6))-(C6/(r6))
430430
elseif r < R
431-
return ((78*C12*(invR6*invR6*invR2)) - (21*C6*(invR2*invR6)))*(r^2) -
431+
return λ * ((78*C12*(invR6*invR6*invR2)) - (21*C6*(invR2*invR6)))*(r^2) -
432432
((168*C12*(invR6*invR6*invR)) - (48*C6*(invR6*invR)))*r +
433433
(91*C12*(invR6*invR6)) - (28*C6*(invR6))
434434
end

‎test/interactions.jl‎

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -47,44 +47,44 @@
4747
inter = LennardJonesSoftCoreBeutler(α=0.3, λ=0.5)
4848
@test isapprox(
4949
force(inter, dr14, a1, a1),
50-
SVector(35.093676538737824, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
50+
SVector(17.546838269368916, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
5151
atol=1e-9u"kJ * mol^-1 * nm^-1",
5252
)
5353
@test isapprox(
5454
force(inter, dr13, a1, a1),
55-
SVector(-1.3236594727846438, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
55+
SVector(-0.6618297363923222, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
5656
atol=1e-9u"kJ * mol^-1 * nm^-1",
5757
)
5858
@test isapprox(
5959
potential_energy(inter, dr14, a1, a1),
60-
29.629058917654785u"kJ * mol^-1";
60+
14.814529458827394u"kJ * mol^-1";
6161
atol=1e-9u"kJ * mol^-1",
6262
)
6363
@test isapprox(
6464
potential_energy(inter, dr13, a1, a1),
65-
-0.11464014233084913u"kJ * mol^-1";
65+
-0.05732007116542457u"kJ * mol^-1";
6666
atol=1e-9u"kJ * mol^-1",
6767
)
6868

6969
inter = LennardJonesSoftCoreGapsys(α=0.85, λ=0.5)
7070
@test isapprox(
7171
force(inter, dr14, a1, a1),
72-
SVector(516.8457758610879, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
72+
SVector(258.42288793054365, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
7373
atol=1e-9u"kJ * mol^-1 * nm^-1",
7474
)
7575
@test isapprox(
7676
force(inter, dr13, a1, a1),
77-
SVector(-1.37550973892212, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
77+
SVector(-0.68775486946106, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
7878
atol=1e-9u"kJ * mol^-1 * nm^-1",
7979
)
8080
@test isapprox(
8181
potential_energy(inter, dr14, a1, a1),
82-
51.814929518108826u"kJ * mol^-1";
82+
45.35853723577515u"kJ * mol^-1";
8383
atol=1e-9u"kJ * mol^-1",
8484
)
8585
@test isapprox(
8686
potential_energy(inter, dr13, a1, a1),
87-
-0.1170417308807374u"kJ * mol^-1";
87+
-0.12971227169036878u"kJ * mol^-1";
8888
atol=1e-9u"kJ * mol^-1",
8989
)
9090

@@ -231,44 +231,44 @@
231231
inter = CoulombSoftCoreBeutler(α=0.3, λ=0.5)
232232
@test isapprox(
233233
force(inter, dr13, a1, a1),
234-
SVector(842.06163558501, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
234+
SVector(421.030817792505, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
235235
atol=1e-5u"kJ * mol^-1 * nm^-1",
236236
)
237237
@test isapprox(
238238
force(inter, dr14, a1, a1),
239-
SVector(57.48819886473348, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
239+
SVector(28.74409943236674, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
240240
atol=1e-5u"kJ * mol^-1 * nm^-1",
241241
)
242242
@test isapprox(
243243
potential_energy(inter, dr13, a1, a1),
244-
345.81678703197457u"kJ * mol^-1";
244+
172.90839351598729u"kJ * mol^-1";
245245
atol=1e-5u"kJ * mol^-1",
246246
)
247247
@test isapprox(
248248
potential_energy(inter, dr14, a1, a1),
249-
634.3822744723304u"kJ * mol^-1";
249+
317.1911372361652u"kJ * mol^-1";
250250
atol=1e-5u"kJ * mol^-1",
251251
)
252252

253253
inter = CoulombSoftCoreGapsys(α=0.3, λ=0.5, σQ=1u"nm")
254254
@test isapprox(
255255
force(inter, dr13, a1, a1),
256-
SVector(731.0108657043696, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
256+
SVector(365.5054328521848, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
257257
atol=1e-5u"kJ * mol^-1 * nm^-1",
258258
)
259259
@test isapprox(
260260
force(inter, dr14, a1, a1),
261-
SVector(1276.8008892849266, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
261+
SVector(638.4004446424633, 0.0, 0.0)u"kJ * mol^-1 * nm^-1";
262262
atol=1e-5u"kJ * mol^-1 * nm^-1",
263263
)
264264
@test isapprox(
265265
potential_energy(inter, dr13, a1, a1),
266-
341.80053925707637u"kJ * mol^-1";
266+
170.90026962853818u"kJ * mol^-1";
267267
atol=1e-5u"kJ * mol^-1",
268268
)
269269
@test isapprox(
270270
potential_energy(inter, dr14, a1, a1),
271-
642.9723025054708u"kJ * mol^-1";
271+
321.4861512527354u"kJ * mol^-1";
272272
atol=1e-5u"kJ * mol^-1",
273273
)
274274

0 commit comments

Comments
 (0)