Merge pull request #10 from dextorious/master

Performance improvements

Author: ChrisRackauckas

src/basic_potentials.jl

@@ -34,7 +34,7 @@ end
 
 function Base.show(stream::IO, pp::GravitationalParameters)
     println(stream, "Gravitational:")
-    print(stream, "\tG:"); show(stream, pp.G); 
+    print(stream, "\tG:"); show(stream, pp.G);
     println(stream)
 end
 
@@ -58,7 +58,7 @@ end
 
 function Base.show(stream::IO, pp::ElectrostaticParameters)
     println(stream, "Electrostatic:")
-    print(stream, "\tk:"); show(stream, pp.k); 
+    print(stream, "\tk:"); show(stream, pp.k);
     println(stream)
 end
 
@@ -72,7 +72,7 @@ end
 
 function Base.show(stream::IO, pp::MagnetostaticParameters)
     println(stream, "Magnetostatic:")
-    print(stream, "\tμ/4π:"); show(stream, pp.μ_4π); 
+    print(stream, "\tμ/4π:"); show(stream, pp.μ_4π);
     println(stream)
 end
 
@@ -93,19 +93,19 @@ function pairwise_lennard_jones_acceleration!(dv,
 
     force = @SVector [0.0, 0.0, 0.0];
     ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
-    
+
     for j ∈ indxs
         if j != i
             rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
-            (rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
-                        
-            if success
+            (rij, r, rij_2) = get_interparticle_distance(ri, rj, pbc)
+
+            if rij_2 < p.R2
                 σ_rij_6 = (p.σ2 / rij_2)^3
                 σ_rij_12 = σ_rij_6^2
-                force += (2 * σ_rij_12 - σ_rij_6 ) * rij / rij_2    
-            end        
+                force += (2 * σ_rij_12 - σ_rij_6 ) * rij / rij_2
+            end
         end
-    end   
+    end
     @. dv +=  24 * p.ϵ * force / ms[i]
 end
 
@@ -119,29 +119,28 @@ function pairwise_electrostatic_acceleration!(dv,
     p::ElectrostaticParameters,
     pbc::BoundaryConditions)
 
-    force = @SVector [0.0, 0.0, 0.0];
+    force = @SVector [0.0, 0.0, 0.0]
     ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
     @inbounds for j = 1:n
         if !in(j, exclude[i])
             rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
-
-            (rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
-            if success
-                force += qs[j] * rij / norm(rij)^3
+            rij, r, r2 = get_interparticle_distance(ri, rj, pbc)
+            if r2 < p.R2
+                force += qs[j] * rij / (r*r2)
             end
         end
-    end    
+    end
     @. dv += p.k * qs[i] * force / ms[i]
 end
 
 
-function gravitational_acceleration!(dv, 
+function gravitational_acceleration!(dv,
     rs,
     i::Integer,
     n::Integer,
     bodies::Vector{<:MassBody},
     p::GravitationalParameters)
-    
+
     accel = @SVector [0.0, 0.0, 0.0];
     ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
     @inbounds for j = 1:n
@@ -151,11 +150,11 @@ function gravitational_acceleration!(dv,
             accel -= p.G * bodies[j].m * rij / norm(rij)^3
         end
     end
-    
+
     @. dv += accel
 end
 
-function magnetostatic_dipdip_acceleration!(dv, 
+function magnetostatic_dipdip_acceleration!(dv,
     rs,
     i::Integer,
     n::Integer,
@@ -177,17 +176,17 @@ function magnetostatic_dipdip_acceleration!(dv,
             force += (mi * mij + mj * mir + r * dot(mi, mj) - 5 * r * mir * mij) / rij4
         end
     end
-    
+
     @. dv += 3 * p.μ_4π * force / bodies[i].m
 end
 
-function harmonic_bond_potential_acceleration!(dv, 
+function harmonic_bond_potential_acceleration!(dv,
     rs,
     i::Integer,
     ms::Vector{<:Real},
     neighbouhoods::Dict{Int,Vector{Tuple{Int,Float64}}},
     p::SPCFwParameters)
-    
+
     force = @SVector [0.0, 0.0, 0.0];
     ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
     @inbounds for (j, k) in neighbouhoods[i]
@@ -197,7 +196,7 @@ function harmonic_bond_potential_acceleration!(dv,
         d = r - p.rOH
         force -= d * k * rij / r
     end
-    
+
     @. dv += force / ms[i]
 end
 
@@ -221,14 +220,14 @@ function valence_angle_potential_acceleration!(dv,
     pc = normalize(cross(rcb, rbaXbc))
 
     cosine = dot(rba, rbc) / (norm(rba) * norm(rbc))
-   
+
     if cosine>1
         cosine = 1
     elseif cosine<-1
         cosine =-1
     end
     aHOH = acos(cosine)
-    
+
     force = - p.ka * (aHOH - p.aHOH)
     force_a = pa * force / norm(rba)
     force_c = pc * force / norm(rbc)
@@ -237,4 +236,4 @@ function valence_angle_potential_acceleration!(dv,
     @. dv[:,a] += force_a / ms[a]
     @. dv[:,b] += force_b / ms[b]
     @. dv[:,c] += force_c / ms[c]
-end
+end

src/boundary_conditions.jl

@@ -11,7 +11,7 @@ Base.start(::PeriodicBoundaryConditions) = 1
 
 Base.done(pbc::PeriodicBoundaryConditions, state) = state > length(pbc.boundary)
 
-function Base.next(pbc::PeriodicBoundaryConditions, state) 
+function Base.next(pbc::PeriodicBoundaryConditions, state)
     pbc.boundary[state], state + 1
 end
 
@@ -31,46 +31,41 @@ struct CubicPeriodicBoundaryConditions{cType <: Real} <: BoundaryConditions
     L::cType
 end
 
-function apply_boundary_conditions!(ri, rj, pbc::PeriodicBoundaryConditions, R2)
-    rij = @MVector [Inf, Inf, Inf]
-    success = false
-    rij2 = 0
-
-    for dx in [0,-pbc[2],pbc[2]], dy in [0,-pbc[4],pbc[4]], dz in [0,-pbc[6],pbc[6]]
-        rij = @MVector [ri[1] - rj[1] + dx, ri[2] - rj[2] + dy, ri[3] - rj[3] + dz]
-        for x in (1, 2, 3) 
-            rij[x] -= (pbc[2x] - pbc[2x - 1]) * div(rij[x], (pbc[2x] - pbc[2x - 1]))
-        end
-        rij2 = dot(rij, rij)
-        if  rij2 < R2
-            success = true
-            break
-        end
-    end
-    return (rij, rij2, success)
+function get_interparticle_distance(ri, rj, pbc::PeriodicBoundaryConditions)
+    rij = ri - rj
+    x, y, z = rij
+    while x  < pbc[1]   x += pbc[2]-pbc[1] end
+    while x >= pbc[2]   x -= pbc[2]-pbc[1] end
+    while y  < pbc[3]   y += pbc[4]-pbc[3] end
+    while y >= pbc[4]   y -= pbc[4]-pbc[3] end
+    while z  < pbc[5]   z += pbc[6]-pbc[5] end
+    while z >= pbc[6]   z -= pbc[6]-pbc[5] end
+    rij = @SVector [x, y, z]
+    r2 = rij[1]^2 + rij[2]^2 + rij[3]^2
+    r = sqrt(r2)
+    return (rij, r, r2)
 end
 
-function apply_boundary_conditions!(ri, rj, pbc::CubicPeriodicBoundaryConditions, R2)
-    rij = @MVector [Inf, Inf, Inf]
-    success = false
-    rij2 = 0
-
-    shifts = @SVector [0,-pbc.L,pbc.L]
-    for dx in shifts, dy in shifts, dz in shifts
-        rij = @MVector [ri[1] - rj[1] + dx, ri[2] - rj[2] + dy, ri[3] - rj[3] + dz]
-        for x in (1, 2, 3) 
-            rij[x] -= pbc.L * div(rij[x], pbc.L)
-        end
-        rij2 = dot(rij, rij)
-        if  rij2 < R2
-            success = true
-            break
-        end
-    end
-    return (rij, rij2, success)
+function get_interparticle_distance(ri, rj, bc::CubicPeriodicBoundaryConditions)
+    rij = ri - rj
+    x, y, z = rij
+    size = bc.L
+    radius = 0.5 * size
+    while x >= radius    x -= size end
+    while x < -radius    x += size end
+    while y >= radius    y -= size end
+    while y < -radius    y += size end
+    while z >= radius    z -= size end
+    while z < -radius    z += size end
+    rij = @SVector [x, y, z]
+    r2 = rij[1]^2 + rij[2]^2 + rij[3]^2
+    r = sqrt(r2)
+    return (rij, r, r2)
 end
 
-function apply_boundary_conditions!(ri, rj, pbc::BoundaryConditions, R2)
+function get_interparticle_distance(ri, rj, ::BoundaryConditions)
     rij = ri - rj
-    (rij, dot(rij,rij),true)
-end
+    r2 = rij[1]^2 + rij[2]^2 + rij[3]^2
+    r = sqrt(r2)
+    (rij, r, r2)
+end

src/nbody_simulation_result.jl

@@ -32,7 +32,7 @@ end
 
 function get_velocity(sr::SimulationResult, time::Real, i::Integer=0)
     n = get_coordinate_vector_count(sr.simulation.system)
-    
+
     if typeof(sr.solution[1]) <: RecursiveArrayTools.ArrayPartition
         velocities = sr(time).x[1]
         if i <= 0
@@ -183,13 +183,13 @@ function lennard_jones_potential(p::LennardJonesParameters, indxs::Vector{<:Inte
             j = indxs[ind_j]
             rj = @SVector [coordinates[1, j], coordinates[2, j], coordinates[3, j]]
 
-            (rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
+            (rij, r, r2) = get_interparticle_distance(ri, rj, pbc)
 
-            if !success
-                rij_2 = p.R2
+            if !(r2 < p.R2)
+                r2 = p.R2
             end
 
-            σ_rij_6 = (p.σ2 / rij_2)^3
+            σ_rij_6 = (p.σ2 / r2)^3
             σ_rij_12 = σ_rij_6^2
             e_lj += (σ_rij_12 - σ_rij_6 )
         end
@@ -199,24 +199,25 @@ function lennard_jones_potential(p::LennardJonesParameters, indxs::Vector{<:Inte
 end
 
 function electrostatic_potential(p::ElectrostaticParameters, indxs::Vector{<:Integer}, exclude::Dict{Int,Vector{Int}}, qs, rs, pbc::BoundaryConditions)
-    e_el = 0
+    e_el = 0.0
 
     n = length(indxs)
     for ind_i = 1:n
         i = indxs[ind_i]
         ri = @SVector [rs[1, i], rs[2, i], rs[3, i]]
-        e_el_i = 0
+        e_el_i = 0.0
         for ind_j = ind_i + 1:n
             j = indxs[ind_j]
             if !in(j, exclude[i])
                 rj = @SVector [rs[1, j], rs[2, j], rs[3, j]]
 
-                (rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, p.R2)
-                if !success
-                    rij = p.R
-                end
+                (rij, r, r2) = get_interparticle_distance(ri, rj, pbc)
 
-                e_el_i += qs[j] / norm(rij)
+                if r2 < p.R2
+                    e_el_i += qs[j] / r
+                else
+                    e_el_i += qs[j] / p.R
+                end
             end
         end
         e_el += e_el_i * qs[i]
@@ -461,11 +462,10 @@ function rdf(sr::SimulationResult)
                 j = indxs[ind_j]
                 rj = @SVector [cc[1, j], cc[2, j], cc[3, j]]
 
-                (rij, rij_2, success) = apply_boundary_conditions!(ri, rj, pbc, (0.5 * pbc.L)^2)
+                (rij, r, r2) = get_interparticle_distance(ri, rj, pbc)
 
-                if success
-                    rij_1 = sqrt(rij_2)
-                    bin = ceil(Int, rij_1 / dr)
+                if r2 < (0.5 * pbc.L)^2
+                    bin = ceil(Int, r / dr)
                     if bin > 1 && bin <= maxbin
                         hist[bin] += 2
                     end
@@ -567,12 +567,12 @@ function write_pdb_data(f::IO, sr::SimulationResult{<:WaterSPCFw})
             @. cc[:, indH2] = cc[:,indO] + cc0[:, indH2] - cc0[:, indO]
         end
 
-        count+=1      
+        count+=1
         println(f,rpad("MODEL",10), count)
         println(f,"REMARK 250 time=$t picoseconds")
         for i ∈ 1:n
             indO, indH1, indH2 = 3 * (i - 1) + 1, 3 * (i - 1) + 2, 3 * (i - 1) + 3
-            
+
             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))
             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))
             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 +588,11 @@ function write_pdb_data(f::IO, sr::SimulationResult)
     for t in sr.solution.t
         cc = 10*get_position(sr, t)
         map!(x ->  x -= L * floor(x / L), cc, cc)
-        count+=1      
+        count+=1
         println(f,rpad("MODEL",10), count)
         println(f,"REMARK 250 time=$t steps")
         for i ∈ 1:n
-            
+
             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))
         end
         println(f,"ENDMDL")
@@ -656,7 +656,7 @@ function extract_from_pdb(file)
             break
         end
     end
-    
+
 
     ln_time = readline(file)
     parts = split(ln_time)
@@ -687,7 +687,7 @@ function extract_from_pdb(file)
 
             push!(wms, WaterMolecule(o,h1,h2))
         end
-    end       
+    end
 
     wms
-end
+end