Merge pull request #98 from m3g/celllistmap-0.10.0

Celllistmap 0.10.0

Author: lmiq

CHANGELOG.md

@@ -12,6 +12,7 @@ MolSimToolkit.jl Changelog
 
 Version 2.0.1-DEV
 --------------
+- ![INFO][badge-info] Requires CellListMap 0.10.0 (and thus PDBTools 3.25.2)
 
 Version 2.0.0
 --------------

Project.toml

@@ -36,7 +36,7 @@ Plotting = "Plots"
 [compat]
 Aqua = "0.8.13"
 BenchmarkTools = "1.3"
-CellListMap = "0.9.11"
+CellListMap = "0.10.0"
 Chemfiles = "0.10.31"
 ChunkSplitters = "3.1.1"
 Dates = "1.10"
@@ -50,7 +50,7 @@ MDLovoFit_jll = "20.1"
 MolSimToolkitShared = "1.5"
 OffsetArrays = "1.13"
 OrderedCollections = "1.8.1"
-PDBTools = "3.21.0"
+PDBTools = "3.25.2"
 Plots = "1.40"
 Printf = "1.10"
 ProgressMeter = "1.10"

src/BlockAverages.jl

@@ -247,7 +247,6 @@ function block_average(
     tau_int = 1 + 2 * sum(_autocor[i] for i in 2:i95-1; init=0.0)
     n_eff = n / tau_int
     xmean_stderr_neff = std(x) / sqrt(n_eff)
-
     return BlockAverageData{TINPUT,DT}(
         x_input,
         xmean * oneunit(TINPUT),

src/MolecularMinimumDistances/AllPairs.jl

@@ -4,6 +4,8 @@ struct AllPairs{T,F1<:Function,F2<:Function} <: SystemPairs
     ymol_indices::F2
 end
 
+_uround(x) = round(x / oneunit(x); digits = 2) * oneunit(x)
+
 import Base.show
 function Base.show(io::IO, ::MIME"text/plain", sys::AllPairs)
     print(io,chomp("""
@@ -42,7 +44,7 @@ more general `xmol_indices` and/or `ymol_indices` functions,
 which, for each atomic index, returns the corresponding
 molecular index (which is `mol_indices(i) = (i-1)%n + 1` where `n` is the
 number of atoms per molecule if all molecules have the same number of
-atoms and are continously stored in the array of positions). 
+atoms and are continuously stored in the array of positions). 
 
 # Examples
 
@@ -105,18 +107,12 @@ function AllPairs(;
 end
 
 #=
-    function update_list!(
-        i, j, d2,
-        mol_indices_i::Fi,
-        mol_indices_j::Fj,
-        lists::Tuple{T,T}
-    ) where {Fi<:Function, Fj<:Function, T<:AbstractVector{<:MinimumDistance}}
-
 
 Function to update the minimum distance in the case where two lists are being constructed.
 
 =#
-function update_list!(i, j, d2, lists, system::AllPairs)
+function update_list!(pair, lists, system::AllPairs)
+    (;i, j, d2) = pair
     x_list = lists[1]
     y_list = lists[2]
     d = sqrt(d2)
@@ -141,11 +137,11 @@ end
         MolSimToolkit.Testing.namd_pdb,
         MolSimToolkit.Testing.namd_traj,
     )
-    first_frame!(simulation)
-    p = positions(current_frame(simulation))
-    uc = unitcell(current_frame(simulation))
-    xsolvent = zeros(eltype(p), length(popc))
-    xsolute = zeros(eltype(p), length(protein))
+    f = first_frame!(simulation)
+    p = positions(f)
+    uc = unitcell(f)
+    xsolvent = p[popc]    
+    xsolute = p[protein]
     sys = AllPairs(
         xpositions = xsolvent,
         ypositions = xsolute,
@@ -159,17 +155,17 @@ end
     for (iframe, frame) in enumerate(simulation)
         pos = positions(frame)
         local uc = unitcell(current_frame(simulation))
-        sys.xpositions .= @view(pos[popc])
-        sys.ypositions .= @view(pos[protein])
-        sys.unitcell = uc.orthorhombic ? diag(uc.matrix) : uc.matrix
+        sys.system.xpositions .= @view(pos[popc])
+        sys.system.ypositions .= @view(pos[protein])
+        sys.system.unitcell = uc.orthorhombic ? diag(uc.matrix) : uc.matrix
         md = minimum_distances!(sys)
         xmd_min[iframe] = minimum(p -> p.d, md[1])
         ymd_indices[iframe] = minimum(p -> p.i, md[2])
         # Test direct (out-of-place) call
         if iframe == 1
             md_out = minimum_distances(
-                xpositions = sys.xpositions,
-                ypositions = sys.ypositions,
+                xpositions = xsolvent,
+                ypositions = xsolute,
                 xn_atoms_per_molecule = 134,
                 yn_atoms_per_molecule = length(protein),
                 cutoff = 6.0,

src/MolecularMinimumDistances/CrossPairs.jl

@@ -38,7 +38,7 @@ more general `mol_indices` function, which, for each atomic index, returns the
 corresponding
 molecular index (which is `mol_indices(i) = (i-1)%n + 1` where `n` is the
 number of atoms per molecule if all molecules have the same number of
-atoms and are continously stored in the array of positions). 
+atoms and are continuously stored in the array of positions). 
 
 # Examples
 
@@ -96,8 +96,8 @@ end
 # Functions that return the atoms of the second set that are closer to
 # each molecule of the first set (only one list is returned).
 #
-function update_list!(i, j, d2, list, system::CrossPairs)
-    d = sqrt(d2)
+function update_list!(pair, list, system::CrossPairs)
+    (; i, j, d) = pair
     imol = system.mol_indices(i)
     if d < list[imol].d
         list[imol] = MinimumDistance(true, i, j, d)
@@ -115,11 +115,11 @@ end
         MolSimToolkit.Testing.namd_pdb,
         MolSimToolkit.Testing.namd_traj,
     )
-    first_frame!(simulation)
-    p = positions(current_frame(simulation))
-    uc = unitcell(current_frame(simulation))
-    xsolvent = zeros(eltype(p), length(popc))
-    xsolute = zeros(eltype(p), length(protein))
+    f = first_frame!(simulation)
+    p = positions(f)
+    uc = unitcell(f)
+    xsolvent = p[popc]
+    xsolute = p[protein]
     sys = CrossPairs(
         xpositions = xsolvent,
         ypositions = xsolute,
@@ -131,19 +131,19 @@ end
     for (iframe, frame) in enumerate(simulation)
         pos = positions(frame)
         local uc = unitcell(frame)
-        sys.xpositions .= @view(pos[popc])
-        sys.ypositions .= @view(pos[protein])
-        sys.unitcell = uc.orthorhombic ? diag(uc.matrix) : uc.matrix
+        sys.system.xpositions .= @view(pos[popc])
+        sys.system.ypositions .= @view(pos[protein])
+        sys.system.unitcell = uc.orthorhombic ? diag(uc.matrix) : uc.matrix
         md = minimum_distances!(sys)
         md_count[iframe] = count(p -> p.within_cutoff, md)
         # Test direct (out-of-place) call
         if iframe == 1
             md_out = minimum_distances(
-                xpositions = sys.xpositions,
-                ypositions = sys.ypositions,
+                xpositions = xsolvent,
+                ypositions = xsolute,
                 xn_atoms_per_molecule = 134,
                 cutoff = 6.0,
-                unitcell = sys.unitcell
+                unitcell = sys.system.unitcell
             )
             @test all(md_out .≈ md)
         end

src/MolecularMinimumDistances/MolecularMinimumDistances.jl

@@ -1,14 +1,12 @@
 module MolecularMinimumDistances
 
-import TestItems: @testitem
-import DocStringExtensions: TYPEDEF, TYPEDFIELDS
+using TestItems: @testitem
+using DocStringExtensions: TYPEDEF, TYPEDFIELDS
 
-import LinearAlgebra: diag
+using LinearAlgebra: diag
+using StaticArrays: SVector
+using CellListMap: ParticleSystem, pairwise!, update!
 import PDBTools: distance
-import StaticArrays: SVector
-import CellListMap: _uround # interal function, used for showing the unitcell
-import CellListMap: ParticleSystem, map_pairwise, map_pairwise!,
-    update_cutoff!, update_unitcell!
 
 export MinimumDistance
 export SelfPairs, CrossPairs, AllPairs
@@ -205,7 +203,7 @@ julia> init_list(x, i -> (i-1)÷2 + 1)
 
 ```
 
-The above annonymous function `i -> (i-1)÷2 + 1` is equivalent to `i -> mol_indices(i,2)`,
+The above anonymous function `i -> (i-1)÷2 + 1` is equivalent to `i -> mol_indices(i,2)`,
 and can be generalized if the the number of atoms of each molecule is not the same.
 
 =#
@@ -231,7 +229,7 @@ init_list(::Type{T}, n::Integer) where {T} = fill(zero(MinimumDistance{T}), n)
 # Note that we have two types of output variables here: List, and a tuple of List.
 # The List is simply an array of `MinimumDistance{T}`, and we have defined above
 # `copy` and `zero` methods for this type, such that we only need to define 
-# that reseting this variable consists of returing its zero, and set up the reducer.
+# that resetting this variable consists of returning its zero, and set up the reducer.
 # The methods for abstract arrays will take care of the rest.
 #
 # For the Tuple of lists, we need to be more explicit, and define appropriate copy_output,
@@ -263,7 +261,7 @@ end
 Function that computes the minimum distances for an initialized system,
 of `SelfPairs`, `CrossPairs`, or `AllPairs` types. 
 
-The function returs a `Vector{MinimumDistance}` cor `SelfPairs` and `CrossPairs`
+The function returns a `Vector{MinimumDistance}` cor `SelfPairs` and `CrossPairs`
 inputs, and a Tuple of two of such vectors for the `AllPairs` input types.
 
 This function is used as an advanced alternative from preallocated system inputs. Only a few allocations 
@@ -303,10 +301,7 @@ julia> @btime minimum_distances!(\$sys);
 
 """
 function minimum_distances!(sys)
-    map_pairwise!(
-        (x, y, i, j, d2, list) -> update_list!(i, j, d2, list, sys),
-        sys.system
-    )
+    pairwise!((pair, list) -> update_list!(pair, list, sys), sys.system)
     return sys.minimum_distances
 end
 
@@ -482,24 +477,12 @@ abstract type SystemPairs end
 
 import Base: getproperty, propertynames
 getproperty(sys::SystemPairs, s::Symbol) = getproperty(sys, Val(s))
+getproperty(sys::SystemPairs, ::Val{:minimum_distances}) = sys.system.output
 getproperty(sys::SystemPairs, ::Val{:system}) = getfield(sys, :system)
 getproperty(sys::SystemPairs, ::Val{:mol_indices}) = getfield(sys, :mol_indices)
 getproperty(sys::SystemPairs, ::Val{:xmol_indices}) = getfield(sys, :xmol_indices)
 getproperty(sys::SystemPairs, ::Val{:ymol_indices}) = getfield(sys, :ymol_indices)
-getproperty(sys::SystemPairs, ::Val{:minimum_distances}) = sys.system.output
-getproperty(sys::SystemPairs, ::Val{:xpositions}) = sys.system.xpositions
-getproperty(sys::SystemPairs, ::Val{:ypositions}) = sys.system.ypositions
-getproperty(sys::SystemPairs, ::Val{:cutoff}) = sys.system.cutoff
-getproperty(sys::SystemPairs, ::Val{:unitcell}) = sys.system.unitcell
-getproperty(sys::SystemPairs, ::Val{:parallel}) = sys.system.parallel
-propertynames(sys::SystemPairs, private::Bool) =
-    (:system, :mol_indices, :minimum_distances, :xpositions, :ypositions, :unitcell, :cutoff)
-
-import Base: setproperty!
-setproperty!(sys::SystemPairs, s::Symbol, value) = setproperty!(sys, Val(s), value)
-setproperty!(sys::SystemPairs, ::Val{:cutoff}, cutoff) = update_cutoff!(sys.system, cutoff)
-setproperty!(sys::SystemPairs, ::Val{:unitcell}, unitcell) = update_unitcell!(sys.system, unitcell)
-setproperty!(sys::SystemPairs, ::Val{:parallel}, parallel) = sys.system.parallel = parallel
+propertynames(sys::SystemPairs, private::Bool) = (:system, :minimum_distances)
 
 #
 # Functions for when the lists of minimum-distances is that of a single

src/MolecularMinimumDistances/SelfPairs.jl

@@ -33,7 +33,7 @@ more general `mol_indices` function, which, for each atomic index, returns the
 corresponding
 molecular index (which is `mol_indices(i) = (i-1)%n + 1` where `n` is the
 number of atoms per molecule if all molecules have the same number of
-atoms and are continously stored in the array of positions). 
+atoms and are continuously stored in the array of positions). 
 
 # Examples
 
@@ -84,10 +84,11 @@ function SelfPairs(;
 end
 
 #
-# This file containst the functions for single-sets, that is for those cases where
+# This file contains the functions for single-sets, that is for those cases where
 # the list of minimum-distances is between the molecules of a single component.
 #
-function update_list!(i, j, d2, list::List, system::SelfPairs)
+function update_list!(pair, list::List, system::SelfPairs)
+    (; i, j, d2) = pair
     imol = system.mol_indices(i)
     jmol = system.mol_indices(j)
     if imol != jmol
@@ -111,10 +112,10 @@ end
         MolSimToolkit.Testing.namd_pdb,
         MolSimToolkit.Testing.namd_traj,
     )
-    first_frame!(simulation)
-    p = positions(current_frame(simulation))
-    uc = unitcell(current_frame(simulation))
-    xsolvent = zeros(eltype(p), length(popc))
+    f = first_frame!(simulation)
+    p = positions(f)
+    uc = unitcell(f)
+    xsolvent = p[popc]
     sys = SelfPairs(
         xpositions = xsolvent,
         cutoff = 6.0,
@@ -125,17 +126,17 @@ end
     for (iframe, frame) in enumerate(simulation)
         pos = positions(frame)
         local uc = unitcell(frame)
-        sys.xpositions .= @view(pos[popc])
-        sys.unitcell = uc.orthorhombic ? diag(uc.matrix) : uc.matrix
+        sys.system.xpositions .= @view(pos[popc])
+        sys.system.unitcell = uc.orthorhombic ? diag(uc.matrix) : uc.matrix
         md = minimum_distances!(sys)
         md_min[iframe] = minimum(p -> p.d, md)
         # Test direct (out-of-place) call
         if iframe == 1
             md_out = minimum_distances(
-                xpositions = sys.xpositions,
+                xpositions = xsolvent,
                 xn_atoms_per_molecule = 134,
                 cutoff = 6.0,
-                unitcell = sys.unitcell
+                unitcell = sys.system.unitcell
             )
             @test all(md_out .≈ md)
         end

src/Reweighting/Reweight.jl

@@ -127,7 +127,7 @@ function reweight(
             output = 0.0,
             output_name = :total_energy
         )
-        energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
+        energy_vec[iframe] = pairwise!((pair, total_energy) -> total_energy + f_perturbation(pair.i, pair.j, pair.d/10), system)
     end
     @. prob_rel_vec = exp(-(energy_vec)/k*T)
     prob_vec = prob_rel_vec/sum(prob_rel_vec)
@@ -159,7 +159,7 @@ function reweight(
             output = 0.0,
             output_name = :total_energy
         )
-        energy_vec[iframe] = map_pairwise!((x, y, i, j, d2, total_energy) -> total_energy + f_perturbation(i, j, sqrt(d2)/10), system)
+        energy_vec[iframe] = pairwise!((pair, total_energy) -> total_energy + f_perturbation(pair.i, pair.j, pair.d/10), system)
     end
     @. prob_rel_vec = exp(-(energy_vec)/k*T)
     prob_vec = prob_rel_vec/sum(prob_rel_vec)

src/Reweighting/Reweighting.jl

@@ -1,7 +1,7 @@
 module Reweighting
 
 using LinearAlgebra: diag
-using CellListMap: ParticleSystem, map_pairwise, map_pairwise!
+using CellListMap: ParticleSystem, pairwise!
 using ..MolSimToolkit: Simulation, positions, unitcell
 
 export reweight, lennard_jones_perturbation, poly_decay_perturbation, gaussian_decay_perturbation

src/hydrogen_bonds/hydrogen_bonds.jl

@@ -1,6 +1,6 @@
 using Base.Threads: @threads
 using ChunkSplitters: index_chunks
-using CellListMap: CellListMap, ParticleSystem, map_pairwise!, update_unitcell!
+using CellListMap: CellListMap, ParticleSystem, pairwise!, update!
 using OrderedCollections: OrderedDict
 
 include("./internals.jl")
@@ -119,7 +119,7 @@ function PDBTools.hydrogen_bonds(
                 if sel1 == sel2
                     s1 = selection_data[sel1]
                     sys.xpositions .= @view(current_positions[s1.inds])
-                    update_unitcell!(sys, uc.matrix)
+                    update!(sys; unitcell=uc.matrix)
                     number_of_hbonds = count_hbonds(
                         sys, s1, angle_cutoff, electronegative_elements
                     )
@@ -128,7 +128,7 @@ function PDBTools.hydrogen_bonds(
                     s2 = selection_data[sel2]
                     sys.xpositions .= @view(current_positions[s1.inds])
                     sys.ypositions .= @view(current_positions[s2.inds])
-                    update_unitcell!(sys, uc.matrix)
+                    update!(sys; unitcell=uc.matrix)
                     number_of_hbonds = count_hbonds(
                         sys, s1, s2, sel1, sel2, angle_cutoff, electronegative_elements
                     )