Skip to content

Commit 85586c9

Browse files
dextoriousdextorious
committed
Reworked remaining BC implementations
1 parent 806281e commit 85586c9

File tree

3 files changed

+76
-76
lines changed

3 files changed

+76
-76
lines changed

‎src/basic_potentials.jl‎

Lines changed: 26 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ end
3434

3535
function Base.show(stream::IO, pp::GravitationalParameters)
3636
println(stream, "Gravitational:")
37-
print(stream, "\tG:"); show(stream, pp.G);
37+
print(stream, "\tG:"); show(stream, pp.G);
3838
println(stream)
3939
end
4040

@@ -58,7 +58,7 @@ end
5858

5959
function Base.show(stream::IO, pp::ElectrostaticParameters)
6060
println(stream, "Electrostatic:")
61-
print(stream, "\tk:"); show(stream, pp.k);
61+
print(stream, "\tk:"); show(stream, pp.k);
6262
println(stream)
6363
end
6464

@@ -72,7 +72,7 @@ end
7272

7373
function Base.show(stream::IO, pp::MagnetostaticParameters)
7474
println(stream, "Magnetostatic:")
75-
print(stream, "\tμ/4π:"); show(stream, pp.μ_4π);
75+
print(stream, "\tμ/4π:"); show(stream, pp.μ_4π);
7676
println(stream)
7777
end
7878

@@ -93,19 +93,19 @@ function pairwise_lennard_jones_acceleration!(dv,
9393

9494
force = @SVector [0.0, 0.0, 0.0];
9595
ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
96-
96+
9797
for j ∈ indxs
9898
if j != i
9999
rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
100-
(rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
101-
102-
if success
100+
(rij, r, rij_2) = get_interparticle_distance(ri, rj, pbc)
101+
102+
if rij_2 < p.R2
103103
σ_rij_6 = (p.σ2 / rij_2)^3
104104
σ_rij_12 = σ_rij_6^2
105-
force += (2 * σ_rij_12 - σ_rij_6 ) * rij / rij_2
106-
end
105+
force += (2 * σ_rij_12 - σ_rij_6 ) * rij / rij_2
106+
end
107107
end
108-
end
108+
end
109109
@. dv += 24 * p.ϵ * force / ms[i]
110110
end
111111

@@ -119,29 +119,28 @@ function pairwise_electrostatic_acceleration!(dv,
119119
p::ElectrostaticParameters,
120120
pbc::BoundaryConditions)
121121

122-
force = @SVector [0.0, 0.0, 0.0];
122+
force = @SVector [0.0, 0.0, 0.0]
123123
ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
124124
@inbounds for j = 1:n
125125
if !in(j, exclude[i])
126126
rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
127-
128-
(rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
129-
if success
130-
force += qs[j] * rij / norm(rij)^3
127+
rij, r, r2 = get_interparticle_distance(ri, rj, pbc)
128+
if r2 < p.R2
129+
force += qs[j] * rij / (r*r2)
131130
end
132131
end
133-
end
132+
end
134133
@. dv += p.k * qs[i] * force / ms[i]
135134
end
136135

137136

138-
function gravitational_acceleration!(dv,
137+
function gravitational_acceleration!(dv,
139138
rs,
140139
i::Integer,
141140
n::Integer,
142141
bodies::Vector{<:MassBody},
143142
p::GravitationalParameters)
144-
143+
145144
accel = @SVector [0.0, 0.0, 0.0];
146145
ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
147146
@inbounds for j = 1:n
@@ -151,11 +150,11 @@ function gravitational_acceleration!(dv,
151150
accel -= p.G * bodies[j].m * rij / norm(rij)^3
152151
end
153152
end
154-
153+
155154
@. dv += accel
156155
end
157156

158-
function magnetostatic_dipdip_acceleration!(dv,
157+
function magnetostatic_dipdip_acceleration!(dv,
159158
rs,
160159
i::Integer,
161160
n::Integer,
@@ -177,17 +176,17 @@ function magnetostatic_dipdip_acceleration!(dv,
177176
force += (mi * mij + mj * mir + r * dot(mi, mj) - 5 * r * mir * mij) / rij4
178177
end
179178
end
180-
179+
181180
@. dv += 3 * p.μ_4π * force / bodies[i].m
182181
end
183182

184-
function harmonic_bond_potential_acceleration!(dv,
183+
function harmonic_bond_potential_acceleration!(dv,
185184
rs,
186185
i::Integer,
187186
ms::Vector{<:Real},
188187
neighbouhoods::Dict{Int,Vector{Tuple{Int,Float64}}},
189188
p::SPCFwParameters)
190-
189+
191190
force = @SVector [0.0, 0.0, 0.0];
192191
ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
193192
@inbounds for (j, k) in neighbouhoods[i]
@@ -197,7 +196,7 @@ function harmonic_bond_potential_acceleration!(dv,
197196
d = r - p.rOH
198197
force -= d * k * rij / r
199198
end
200-
199+
201200
@. dv += force / ms[i]
202201
end
203202

@@ -221,14 +220,14 @@ function valence_angle_potential_acceleration!(dv,
221220
pc = normalize(cross(rcb, rbaXbc))
222221

223222
cosine = dot(rba, rbc) / (norm(rba) * norm(rbc))
224-
223+
225224
if cosine>1
226225
cosine = 1
227226
elseif cosine<-1
228227
cosine =-1
229228
end
230229
aHOH = acos(cosine)
231-
230+
232231
force = - p.ka * (aHOH - p.aHOH)
233232
force_a = pa * force / norm(rba)
234233
force_c = pc * force / norm(rbc)
@@ -237,4 +236,4 @@ function valence_angle_potential_acceleration!(dv,
237236
@. dv[:,a] += force_a / ms[a]
238237
@. dv[:,b] += force_b / ms[b]
239238
@. dv[:,c] += force_c / ms[c]
240-
end
239+
end

‎src/boundary_conditions.jl‎

Lines changed: 32 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -31,40 +31,42 @@ struct CubicPeriodicBoundaryConditions{cType <: Real} <: BoundaryConditions
3131
L::cType
3232
end
3333

34-
function apply_boundary_conditions!(ri, rj, pbc::PeriodicBoundaryConditions, R2)
35-
rij = @MVector [Inf, Inf, Inf]
36-
success = false
37-
rij2 = 0
38-
39-
for dx in [0,-pbc[2],pbc[2]], dy in [0,-pbc[4],pbc[4]], dz in [0,-pbc[6],pbc[6]]
40-
rij = @MVector [ri[1] - rj[1] + dx, ri[2] - rj[2] + dy, ri[3] - rj[3] + dz]
41-
for x in (1, 2, 3)
42-
rij[x] -= (pbc[2x] - pbc[2x - 1]) * div(rij[x], (pbc[2x] - pbc[2x - 1]))
43-
end
44-
rij2 = dot(rij, rij)
45-
if rij2 < R2
46-
success = true
47-
break
48-
end
49-
end
50-
return (rij, rij2, success)
34+
function get_interparticle_distance(ri, rj, pbc::PeriodicBoundaryConditions)
35+
rij = ri - rj
36+
r, r2 = zero(eltype(ri)), zero(eltype(ri))
37+
x, y, z = rij
38+
while x < pbc[1] x += pbc[2]-pbc[1] end
39+
while x >= pbc[2] x -= pbc[2]-pbc[1] end
40+
while y < pbc[3] y += pbc[4]-pbc[3] end
41+
while y >= pbc[4] y -= pbc[4]-pbc[3] end
42+
while z < pbc[5] z += pbc[6]-pbc[5] end
43+
while z >= pbc[6] z -= pbc[6]-pbc[5] end
44+
rij = @SVector [x, y, z]
45+
r2 = dot(rij, rij)
46+
r = sqrt(r2)
47+
return (rij, r, r2)
5148
end
5249

53-
function apply_boundary_conditions!(ri, rj, pbc::CubicPeriodicBoundaryConditions, R2)
50+
function get_interparticle_distance(ri, rj, bc::CubicPeriodicBoundaryConditions)
5451
rij = ri - rj
55-
x, y, z = rij[1], rij[2], rij[3]
56-
while x >= pbc.L x -= pbc.L end
57-
while x < -pbc.L x += pbc.L end
58-
while y >= pbc.L y -= pbc.L end
59-
while y < -pbc.L y += pbc.L end
60-
while z >= pbc.L z -= pbc.L end
61-
while z < -pbc.L z += pbc.L end
62-
rij = SVector{3,eltype(R2)}(x, y, z)
63-
rij2 = dot(rij, rij)
64-
return (rij, rij2, rij2 < R2)
52+
x, y, z = rij
53+
size = bc.L
54+
radius = 0.5 * size
55+
while x >= radius x -= size end
56+
while x < -radius x += size end
57+
while y >= radius y -= size end
58+
while y < -radius y += size end
59+
while z >= radius z -= size end
60+
while z < -radius z += size end
61+
rij = @SVector [x, y, z]
62+
r2 = dot(rij, rij)
63+
r = sqrt(r2)
64+
return (rij, r, r2)
6565
end
6666

67-
function apply_boundary_conditions!(ri, rj, pbc::BoundaryConditions, R2)
67+
function get_interparticle_distance(ri, rj, ::BoundaryConditions)
6868
rij = ri - rj
69-
(rij, dot(rij,rij),true)
69+
r2 = dot(rij, rij)
70+
r = sqrt(r2)
71+
(rij, r, r2)
7072
end

‎src/nbody_simulation_result.jl‎

Lines changed: 18 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ end
3232

3333
function get_velocity(sr::SimulationResult, time::Real, i::Integer=0)
3434
n = get_coordinate_vector_count(sr.simulation.system)
35-
35+
3636
if typeof(sr.solution[1]) <: RecursiveArrayTools.ArrayPartition
3737
velocities = sr(time).x[1]
3838
if i <= 0
@@ -183,13 +183,13 @@ function lennard_jones_potential(p::LennardJonesParameters, indxs::Vector{<:Inte
183183
j = indxs[ind_j]
184184
rj = @SVector [coordinates[1, j], coordinates[2, j], coordinates[3, j]]
185185

186-
(rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
186+
(rij, r, r2) = get_interparticle_distance(ri, rj, pbc)
187187

188-
if !success
189-
rij_2 = p.R2
188+
if !(r2 < p.R2)
189+
r2 = p.R2
190190
end
191191

192-
σ_rij_6 = (p.σ2 / rij_2)^3
192+
σ_rij_6 = (p.σ2 / r2)^3
193193
σ_rij_12 = σ_rij_6^2
194194
e_lj += (σ_rij_12 - σ_rij_6 )
195195
end
@@ -211,9 +211,9 @@ function electrostatic_potential(p::ElectrostaticParameters, indxs::Vector{<:Int
211211
if !in(j, exclude[i])
212212
rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
213213

214-
(rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
215-
if !success
216-
rij = p.R
214+
(rij, r, r2) = get_interparticle_distance(ri, rj, pbc)
215+
if !(r2 < p.R2)
216+
rij = p.R # NOTE: hmm?
217217
end
218218

219219
e_el_i += qs[j] / norm(rij)
@@ -461,11 +461,10 @@ function rdf(sr::SimulationResult)
461461
j = indxs[ind_j]
462462
rj = @SVector [cc[1, j], cc[2, j], cc[3, j]]
463463

464-
(rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, (0.5 * pbc.L)^2)
464+
(rij, r, r2) = get_interparticle_distance(ri, rj, pbc)
465465

466-
if success
467-
rij_1 = sqrt(rij_2)
468-
bin = ceil(Int, rij_1 / dr)
466+
if r2 < (0.5 * pbc.L)^2
467+
bin = ceil(Int, r / dr)
469468
if bin > 1 && bin <= maxbin
470469
hist[bin] += 2
471470
end
@@ -567,12 +566,12 @@ function write_pdb_data(f::IO, sr::SimulationResult{<:WaterSPCFw})
567566
@. cc[:, indH2] = cc[:,indO] + cc0[:, indH2] - cc0[:, indO]
568567
end
569568

570-
count+=1
569+
count+=1
571570
println(f,rpad("MODEL",10), count)
572571
println(f,"REMARK 250 time=$t picoseconds")
573572
for i ∈ 1:n
574573
indO, indH1, indH2 = 3 * (i - 1) + 1, 3 * (i - 1) + 2, 3 * (i - 1) + 3
575-
574+
576575
println(f,"HETATM",lpad(indO,5)," ",rpad("O",4),"HOH",lpad(i,6)," ", lpad(@sprintf("%8.3f",cc[1,indO]),8), lpad(@sprintf("%8.3f",cc[2,indO]),8), lpad(@sprintf("%8.3f",cc[3,indO]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("O",2))
577576
println(f,"HETATM",lpad(indH1,5)," ",rpad("H1",4),"HOH",lpad(i,6)," ", lpad(@sprintf("%8.3f",cc[1,indH1]),8), lpad(@sprintf("%8.3f",cc[2,indH1]),8), lpad(@sprintf("%8.3f",cc[3,indH1]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("H",2))
578577
println(f,"HETATM",lpad(indH2,5)," ",rpad("H2",4),"HOH",lpad(i,6)," ", lpad(@sprintf("%8.3f",cc[1,indH2]),8), lpad(@sprintf("%8.3f",cc[2,indH2]),8), lpad(@sprintf("%8.3f",cc[3,indH2]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("H",2))
@@ -588,11 +587,11 @@ function write_pdb_data(f::IO, sr::SimulationResult)
588587
for t in sr.solution.t
589588
cc = 10*get_position(sr, t)
590589
map!(x -> x -= L * floor(x / L), cc, cc)
591-
count+=1
590+
count+=1
592591
println(f,rpad("MODEL",10), count)
593592
println(f,"REMARK 250 time=$t steps")
594593
for i ∈ 1:n
595-
594+
596595
println(f,"HETATM",lpad(i,5)," ",rpad("Ar",4),"Ar",lpad(i,6)," ", lpad(@sprintf("%8.3f",cc[1,i]),8), lpad(@sprintf("%8.3f",cc[2,i]),8), lpad(@sprintf("%8.3f",cc[3,i]),8),lpad("1.00",6),lpad("0.00",6), lpad("",10), lpad("Ar",2))
597596
end
598597
println(f,"ENDMDL")
@@ -656,7 +655,7 @@ function extract_from_pdb(file)
656655
break
657656
end
658657
end
659-
658+
660659

661660
ln_time = readline(file)
662661
parts = split(ln_time)
@@ -687,7 +686,7 @@ function extract_from_pdb(file)
687686

688687
push!(wms, WaterMolecule(o,h1,h2))
689688
end
690-
end
689+
end
691690

692691
wms
693-
end
692+
end

0 commit comments

Comments
 (0)