Skip to content

Commit

Permalink
support nothing as unitcell for non-periodic systems
Browse files Browse the repository at this point in the history
  • Loading branch information
@lmiq
lmiq committed May 1, 2024
1 parent 0d5cea0 commit 499df3f
Show file tree
Hide file tree
Showing 3 changed files with 224 additions and 98 deletions.
2 changes: 1 addition & 1 deletion src/CellListMap.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,8 @@ include("./linearalgebra.jl")
include("./show.jl")
include("./Box.jl")
include("./CellLists.jl")
include("./PeriodicSystems.jl")
include("./CellOperations.jl")
include("./PeriodicSystems.jl")
include("./CoreComputing.jl")

"""
Expand Down
114 changes: 90 additions & 24 deletions src/PeriodicSystems.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,10 @@ using StaticArrays: SVector, FieldVector

import ..CellListMap
using ..CellListMap: INTERNAL
using ..CellListMap: Box, update_box
using ..CellListMap: Box, update_box, limits
using ..CellListMap: CellListPair, Swapped, NotSwapped, unitcelltype
using ..CellListMap: argon_pdb_file # for doc tests
using ..CellListMap: OrthorhombicCell, TriclinicCell, NonPeriodicCell

export PeriodicSystem
export map_pairwise!, map_pairwise
Expand All @@ -29,13 +30,13 @@ const SupportedTypes = Union{Number,SVector,FieldVector}
const SupportedCoordinatesTypes = Union{Nothing,AbstractVector{<:AbstractVector},AbstractMatrix}

"""
PeriodicSystem(
PeriodicSystem(;
xpositions::Union{AbstractVector{<:AbstractVector},AbstractMatrix},
#or
xpositions::Union{AbstractVector{<:AbstractVector},AbstractMatrix},
ypositions::Union{AbstractVector{<:AbstractVector},AbstractMatrix},
# and
unitcell::AbstractVecOrMat,
unitcell::Union{Nothing,AbstractVecOrMat} = nothing,
cutoff::Number,
output::Any;
output_name::Symbol,
Expand All @@ -61,6 +62,9 @@ Constructor of the `PeriodicSystem` type given the positions of the particles.
The unit cell (either a vector for `Orthorhombic` cells or a
full unit cell matrix for `Triclinic` cells), the cutoff used for the
construction of the cell lists and the output variable of the calculations.
If unitcell == nothing, the system is considered not-periodic, in which
case artificial periodic boundaries will be built such that images
are farther from each other than the cutoff.
`output_name` can be set to a symbol that best identifies the output variable.
For instance, if `output_name=:forces`, the forces can be retrieved from the
Expand Down Expand Up @@ -125,7 +129,7 @@ function PeriodicSystem(;
positions::SupportedCoordinatesTypes=nothing,
xpositions::SupportedCoordinatesTypes=nothing,
ypositions::SupportedCoordinatesTypes=nothing,
unitcell::AbstractVecOrMat,
unitcell::Union{AbstractVecOrMat,Nothing}=nothing,
cutoff::Number,
output::Any,
output_name::Symbol=:output,
Expand All @@ -149,22 +153,25 @@ function PeriodicSystem(;
end
input_array = reinterpret(reshape, SVector{dim,eltype(input_array)}, input_array)
end
DIM = if eltype(input_array) isa SVector
length(eltype(input_array))
else
if length(input_array) == 0
# If the array is empty, we cannot determine the dimension, so we assume it is the same as the unit cell
size(unitcell, 1)
if !isnothing(unitcell)
DIM = if eltype(input_array) isa SVector
length(eltype(input_array))
else
length(input_array[1])
if length(input_array) == 0
# If the array is empty, we cannot determine the dimension, so we assume it is the same as the unit cell
size(unitcell, 1)
else
length(input_array[1])
end
end
if DIM != size(unitcell, 1)
throw(DimensionMismatch("Dimension of the unit cell ($(size(unitcell, 1))) must match the dimension of the coordinates ($(length(eltype(input_array))))"))
end
end
if DIM != size(unitcell, 1)
throw(DimensionMismatch("Dimension of the unit cell ($(size(unitcell, 1))) must match the dimension of the coordinates ($(length(eltype(input_array))))"))
end
end
# Single set of positions
if isnothing(ypositions)
unitcell = isnothing(unitcell) ? limits(xpositions) : unitcell
_box = CellListMap.Box(unitcell, cutoff, lcell=lcell)
_cell_list = CellListMap.CellList(xpositions, _box; parallel=parallel, nbatches=nbatches)
_aux = CellListMap.AuxThreaded(_cell_list)
Expand All @@ -173,6 +180,7 @@ function PeriodicSystem(;
sys = PeriodicSystem1{output_name}(xpositions, output, _box, _cell_list, _output_threaded, _aux, parallel)
# Two sets of positions
else
unitcell = isnothing(unitcell) ? limits(xpositions, ypositions) : unitcell
_box = CellListMap.Box(unitcell, cutoff, lcell=lcell)
_cell_list = CellListMap.CellList(xpositions, ypositions, _box; parallel=parallel, nbatches=nbatches, autoswap=autoswap)
_aux = CellListMap.AuxThreaded(_cell_list)
Expand Down Expand Up @@ -265,6 +273,7 @@ CellListMap.unitcelltype(sys::AbstractPeriodicSystem) = unitcelltype(sys._box)
x = rand(SVector{3,Float64}, 100)
@test unitcelltype(PeriodicSystem(positions=x, cutoff=0.1, unitcell=[1, 1, 1], output=0.0)) == OrthorhombicCell
@test unitcelltype(PeriodicSystem(positions=x, cutoff=0.1, unitcell=[1 0 0; 0 1 0; 0 0 1], output=0.0)) == TriclinicCell
@test unitcelltype(PeriodicSystem(positions=x, cutoff=0.1, output=0.0)) == NonPeriodicCell

# Argument errors
@test_throws ArgumentError PeriodicSystem(
Expand Down Expand Up @@ -519,14 +528,14 @@ reset_output!(x::T) where {T<:SupportedTypes} = zero(x)
reset_output!(x::AbstractVecOrMat{T}) where {T} = fill!(x, reset_output!(x[begin]))
const reset_output = reset_output!

"""
#=
_reset_all_output!(output, output_threaded)
$(INTERNAL)
Function that resets the output variable and the threaded copies of it.
"""
=#
function _reset_all_output!(output, output_threaded)
output = reset_output!(output)
for i in eachindex(output_threaded)
Expand Down Expand Up @@ -708,7 +717,7 @@ end
update_unitcell!(system, unitcell)
Function to update the unit cell of the system. The `unicell` must be of the
same type (`OrthorhombicCell` or `TriclinicCell`) of the original `system`
same type (`OrthorhombicCell`, `TriclinicCell`) of the original `system`
(changing the type of unit cell requires reconstructing the system).
The `unitcell` can be a `N×N` matrix or a vector of dimension `N`, where
Expand All @@ -717,9 +726,12 @@ The `unitcell` can be a `N×N` matrix or a vector of dimension `N`, where
This function can be used to update the system geometry in iterative schemes,
where the size of the simulation box changes during the simulation.
!!! note
Manual updating of the unit cell of non-periodic systems is not allowed.
# Example
```jldoctest; filter = r"batches.*" => ""
```jldoctest filter = r"batches.*" => ""
julia> using CellListMap.PeriodicSystems, StaticArrays, PDBTools
julia> xpositions = coor(readPDB(PeriodicSystems.argon_pdb_file));
Expand Down Expand Up @@ -752,6 +764,11 @@ PeriodicSystem1{output} of dimension 3, composed of:
"""
function update_unitcell!(sys, unitcell)
if unitcelltype(sys) == NonPeriodicCell
throw(ArgumentError("""\n
Manual updating of the unit cell of non-periodic systems is not allowed.
"""))
end
sys._box = update_box(sys._box; unitcell=unitcell)
return sys
end
Expand All @@ -773,6 +790,11 @@ end
@test diag(sys2.unitcell) == [2, 2, 2]
a = @ballocated update_unitcell!($sys2, SVector(2, 2, 2)) evals = 1 samples = 1
@test a == 0
# Test throwing error on updating non-periodic unit cells
sys = PeriodicSystem(xpositions=x, cutoff=0.1, output=0.0)
@test_throws ArgumentError update_unitcell!(sys, [1, 1, 1])
sys = PeriodicSystem(xpositions=x, ypositions=y, cutoff=0.1, output=0.0)
@test_throws ArgumentError update_unitcell!(sys, [1, 1, 1])
end

"""
Expand Down Expand Up @@ -812,7 +834,17 @@ PeriodicSystem1 of dimension 3, composed of:
Type of output variable: Float64
```
"""
function update_cutoff!(sys, cutoff)
function update_cutoff!(sys::PeriodicSystem1, cutoff)
if unitcelltype(sys) == NonPeriodicCell
sys._box = Box(limits(sys.xpositions), cutoff)
end
sys._box = update_box(sys._box; cutoff=cutoff)
return sys
end
function update_cutoff!(sys::PeriodicSystem2, cutoff)
if unitcelltype(sys) == NonPeriodicCell
sys._box = Box(limits(sys.xpositions, sys.ypositions), cutoff)
end
sys._box = update_box(sys._box; cutoff=cutoff)
return sys
end
Expand All @@ -821,6 +853,7 @@ end
using BenchmarkTools
using StaticArrays
using CellListMap.PeriodicSystems
using PDBTools
x = rand(SVector{3,Float64}, 1000)
sys1 = PeriodicSystem(xpositions=x, unitcell=[1, 1, 1], cutoff=0.1, output=0.0)
update_cutoff!(sys1, 0.2)
Expand All @@ -833,20 +866,37 @@ end
@test sys2.cutoff == 0.2
a = @ballocated update_cutoff!($sys2, 0.1) evals = 1 samples = 1
@test a == 0

# Update cutoff of non-periodic systems
x = coor(readPDB(PeriodicSystems.argon_pdb_file))
sys1 = PeriodicSystem(xpositions=x, cutoff=8.0, output=0.0)
@test unitcelltype(sys1) == NonPeriodicCell
@test sys1.unitcell ≈ [ 35.63 0.0 0.0; 0.0 35.76 0.0; 0.0 0.0 35.79 ] atol = 1e-2
update_cutoff!(sys1, 10.0)
@test sys1.unitcell ≈ [ 39.83 0.0 0.0; 0.0 39.96 0.0; 0.0 0.0 39.99 ] atol = 1e-2
a = @ballocated update_cutoff!($sys1, 8.0) evals = 1 samples = 1
@test a == 0
sys2 = PeriodicSystem(xpositions=x[1:50], ypositions=x[51:100], cutoff=8.0, output=0.0)
@test unitcelltype(sys2) == NonPeriodicCell
@test sys2.unitcell ≈ [ 35.63 0.0 0.0; 0.0 35.76 0.0; 0.0 0.0 35.79 ] atol = 1e-2
update_cutoff!(sys2, 10.0)
@test sys2.unitcell ≈ [ 39.83 0.0 0.0; 0.0 39.96 0.0; 0.0 0.0 39.99 ] atol = 1e-2
a = @ballocated update_cutoff!($sys2, 8.0) evals = 1 samples = 1
@test a == 0
end

#
# This interface is needed to generate random particle coordinates in ComplexMixtures.jl
#
"""
#=
get_computing_box(sys::AbstractPeriodicSystem)
$(INTERNAL)
Retrieves the computing box of the system. The computing box is large enough to
contain all coordinates of the particles, plus the cutoff.
"""
=#
get_computing_box(sys::AbstractPeriodicSystem) = sys._box.computing_box
@testitem "get_computing_box" begin
using StaticArrays
Expand All @@ -856,16 +906,19 @@ get_computing_box(sys::AbstractPeriodicSystem) = sys._box.computing_box
@test PeriodicSystems.get_computing_box(sys) == ([-0.1, -0.1, -0.1], [1.1, 1.1, 1.1])
end

"""
#=
UpdatePeriodicSystem!
$(INTERNAL)
Updates the cell lists for periodic systems.
"""
=#
function UpdatePeriodicSystem!(sys::PeriodicSystem1, update_lists::Bool=true)
if update_lists
if unitcelltype(sys) == NonPeriodicCell
sys._box = Box(limits(sys.xpositions), sys.cutoff)
end
sys._cell_list = CellListMap.UpdateCellList!(
sys.xpositions,
sys._box,
Expand All @@ -881,11 +934,14 @@ function _update_ref_positions!(cl::CellListPair{V,N,T,Swap}, sys) where {V,N,T,
resize!(cl.ref, length(sys.xpositions))
cl.ref .= sys.xpositions
end
function _update_ref_positions!(cl::CellListPair{V,N,T,Swap}, sys) where {V,N,T,Swap<:Swapped}
function _update_ref_positions!(::CellListPair{V,N,T,Swap}, sys) where {V,N,T,Swap<:Swapped}
throw(ArgumentError("update_lists === false requires autoswap == false for 2-set systems."))
end
function UpdatePeriodicSystem!(sys::PeriodicSystem2, update_lists::Bool=true)
if update_lists
if unitcelltype(sys) == NonPeriodicCell
sys._box = Box(limits(sys.xpositions, sys.ypositions), sys.cutoff)
end
sys._cell_list = CellListMap.UpdateCellList!(
sys.xpositions,
sys.ypositions,
Expand Down Expand Up @@ -931,6 +987,16 @@ end
a = @ballocated PeriodicSystems.UpdatePeriodicSystem!($sys) samples = 1 evals = 1
@test a == 0

# Update non-periodic system
x = rand(SVector{3,Float64}, 1000)
sys = PeriodicSystem(xpositions=x, cutoff=0.1, output=0.0, parallel=false)
a = @ballocated PeriodicSystems.UpdatePeriodicSystem!($sys) samples = 1 evals = 1
@test a == 0
y = rand(SVector{3,Float64}, 1000)
sys = PeriodicSystem(xpositions=x, ypositions=y, cutoff=0.1, output=0.0, parallel=false)
a = @ballocated PeriodicSystems.UpdatePeriodicSystem!($sys) samples = 1 evals = 1
@test a == 0

end

"""
Expand Down
Loading

0 comments on commit 499df3f

Please sign in to comment.