add set_reference_coordinates! function

Project.toml

@@ -15,11 +15,12 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 
 [compat]
-Aqua = "0.6, 0.7, 0.8"
-CellListMap = "0.9"
-LinearAlgebra = "1.9"
+Aqua = "0.8"
+BenchmarkTools = "1.5"
+CellListMap = "0.9.6"
 FiniteDifferences = "0.12.30"
 ForwardDiff = "0.10"
+LinearAlgebra = "1.9"
 NativeFileDialog = "0.2.1"
 PDBTools = "1.5.1"
 Packmol_jll = "20"
@@ -28,12 +29,14 @@ SPGBox = "0.5, 0.6, 0.7"
 StaticArrays = "1"
 Statistics = "1.9"
 Test = "1.9"
-TestItems = "0.1, 1"
 TestItemRunner = "1"
+TestItems = "1"
 julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -42,4 +45,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Aqua", "Test", "Random", "StaticArrays", "ForwardDiff", "FiniteDifferences", "TestItemRunner"]
+test = ["Aqua", "Test", "Random", "StaticArrays", "ForwardDiff", "FiniteDifferences", "TestItemRunner", "CellListMap", "BenchmarkTools"]

src/Packmol.jl

@@ -1,12 +1,13 @@
 module Packmol
 
-using TestItems: @testitem
-using StaticArrays: SVector, SMatrix, @SMatrix
-using LinearAlgebra: norm
+using TestItems: @testitem, @testsnippet
+using StaticArrays: SVector, SMatrix, @SMatrix, MMatrix
+using LinearAlgebra: norm, eigen
 using Statistics: mean
 using Base: @kwdef
 using Base.Threads: @spawn
 using PDBTools: Atom, readPDB, coor 
+import CellListMap
 
 const src_dir = @__DIR__
 

src/rigid_body/rigid_body.jl

@@ -68,8 +68,7 @@ function move!(x::AbstractVector{T}, newcm::T, beta, gamma, theta) where {T<:SVe
     return x
 end
 
-@testitem "move!" begin
-    using StaticArrays
+@testitem "move!" setup=[RigidMolecules] begin
     x = [SVector(1.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0)]
     @test Packmol.move!(copy(x), SVector(0.0, 0.0, 0.0), 0.0, 0.0, 0.0) ≈
           SVector{3,Float64}[[0.5, 0.0, 0.0], [-0.5, 0.0, 0.0]]
@@ -106,49 +105,25 @@ function random_move!(
     # To avoid boundary problems, the center of coordinates are generated in a 
     # much larger region, and wrapped aftwerwards
     scale = 100.0
-
     # Generate random coordinates for the center of mass
     cmin, cmax = CellListMap.get_computing_box(system)
     newcm = SVector{3}(scale * (cmin[i] + rand(RNG, Float64) * (cmax[i] - cmin[i])) for i in 1:3)
-
     # Generate random rotation angles 
     beta = 2π * rand(RNG, Float64)
     gamma = 2π * rand(RNG, Float64)
     theta = 2π * rand(RNG, Float64)
-
     # Take care that this molecule is not split by periodic boundary conditions, by
     # wrapping its coordinates around its reference atom
     for iat in eachindex(x)
         x[iat] = CellListMap.wrap_relative_to(x[iat], x[irefatom], system.unitcell)
     end
-
     # Move molecule to new position
     move!(x, newcm, beta, gamma, theta)
-
     return x
 end
 
-@testitem "random_move!" begin
-    using Packmol
-    using StaticArrays
-    using LinearAlgebra: norm
-    using CellListMap
-    import Random
-
-    function check_internal_distances(x, y)
-        for i = firstindex(x):lastindex(x)-1
-            for j = i+1:lastindex(x)
-                x_ij = norm(x[i] - x[j])
-                y_ij = norm(y[i] - y[j])
-                if !isapprox(x_ij, y_ij)
-                    return false, x_ij, y_ij
-                end
-            end
-        end
-        return true
-    end
 
-    RNG = Random.Xoshiro()
+@testitem "random_move!" setup=[RigidMolecules] begin
     # Orthorhombic cell
     x = [-1.0 .+ 2 * rand(SVector{3,Float64}) for _ = 1:5]
     system = ParticleSystem(positions=x, cutoff=0.1, unitcell=SVector(10.0, 10.0, 10.0), output=0.0)
@@ -157,7 +132,6 @@ end
     @test check_internal_distances(x, Packmol.random_move!(copy(x), 1, system, RNG))
     system.xpositions .= [4.0 .+ 2 * rand(SVector{3,Float64}) for _ = 1:5]
     @test check_internal_distances(x, Packmol.random_move!(copy(x), 1, system, RNG))
-
     # Triclinic cell
     x = [-1.0 .+ 2 * rand(SVector{3,Float64}) for _ = 1:5]
     system = ParticleSystem(positions=x, cutoff=0.1, unitcell=@SMatrix[10.0 5.0 0.0; 0.0 10.0 0.0; 0.0 0.0 10.0], output=0.0)
@@ -166,5 +140,83 @@ end
     @test check_internal_distances(x, Packmol.random_move!(copy(x), 1, system, RNG))
     system.xpositions .= [4.0 .+ 2 * rand(SVector{3,Float64}) for _ = 1:5]
     @test check_internal_distances(x, Packmol.random_move!(copy(x), 1, system, RNG))
+end
+
+#
+# The following function takes the vector of coordinates of a molecule, and rotates
+# and translates it to put the center of mass in the origin and the principal moment
+# of inertia (assuming all equal masses) along the z-axis. That will define the orientation
+# of the reference coordinates.
+#
+function set_reference_coordinates!(x::AbstractVector{<:SVector{3,T}}) where {T<:Real}
+    # Center of mass
+    cm = mean(x)
+    x .= x .- Ref(cm)
+    # Inertia tensor
+    I = zeros(MMatrix{3,3,T,9})
+    for xi in x
+        I[1,1] += xi[2]^2 + xi[3]^2
+        I[2,2] += xi[1]^2 + xi[3]^2
+        I[3,3] += xi[1]^2 + xi[2]^2
+        I[1,2] -= xi[1] * xi[2]
+        I[1,3] -= xi[1] * xi[3]
+        I[2,3] -= xi[2] * xi[3]
+    end
+    I[2,1] = I[1,2]
+    I[3,1] = I[1,3]
+    I[3,2] = I[2,3]
+    I = SMatrix(I)
+    # Diagonalize
+    evals, evecs = eigen(I)
+    # Sort eigenvectors by eigenvalues
+    idx = sortperm(evals, rev=true)
+    evecs = evecs[:, idx]
+    # Rotate
+    for i in eachindex(x)
+        x[i] = evecs' * x[i]
+    end
+    return x
+end
+
+@testitem "set_reference_coordinates!" setup=[rigid_molecules] begin
+    using PDBTools
+    using StaticArrays
+    using BenchmarkTools
+    x = [ SVector(0.0, 0.0, 0.0), SVector(√3/3, √3/3, √3/3) ]
+    Packmol.set_reference_coordinates!(x)
+    @test x[1] ≈ SVector(0., 0., -0.5)
+    @test x[2] ≈ SVector(0., 0., 0.5)
+    x = [ SVector(√3/3, √3/3, √3/3), SVector(0.0, 0.0, 0.0) ]
+    Packmol.set_reference_coordinates!(x)
+    @test x[1] ≈ SVector(0.0, 0.0, 0.5)
+    @test x[2] ≈ SVector(0.0, 0.0, -0.5)
+    water = coor(readPDB(joinpath(Packmol.src_dir, "..", "test", "data", "water.pdb")))
+    water_save = copy(water)
+    a = @ballocated Packmol.set_reference_coordinates!($water) samples=1 evals=1
+    @test a == 0
+    @test mean(water) ≈ SVector(0.0, 0.0, 0.0) atol = 1e-10
+    @test all(isapprox(v[1],0.0,atol=1e-10) for v in water)
+    @test check_internal_distances(water, water_save)
+end
 
-end
+@testsnippet RigidMolecules begin
+    using Packmol
+    using CellListMap
+    using StaticArrays
+    using LinearAlgebra: norm
+    import Random
+    RNG = Random.Xoshiro()
+    # Function that checks that the internal distances of a molecule are preserved
+    function check_internal_distances(x, y)
+        for i = firstindex(x):lastindex(x)-1
+            for j = i+1:lastindex(x)
+                x_ij = norm(x[i] - x[j])
+                y_ij = norm(y[i] - y[j])
+                if !isapprox(x_ij, y_ij)
+                    return false, x_ij, y_ij
+                end
+            end
+        end
+        return true
+    end
+end