move gradient to fg structure

src/data_structures/PackmolSystem.jl

@@ -50,7 +50,6 @@ Base.zero(::Type{MoleculePosition{N,T}}) where {N,T} = MoleculePosition(zero(SVe
     chkgrad::Bool = false
     atoms::Vector{AtomData{D,T}} = AtomData{D,T}[]
     molecule_positions::MoleculePosition{D,T} = MoleculePosition{D,T}[]
-    gradient::MoleculePositions{D,T} = MoleculePosition{D,T}[]
 end
 
 function _indent(s::AbstractString; n=4)

src/interatomic_distance_fg.jl

@@ -5,30 +5,36 @@ using SPGBox: spgbox!
 # in a single pass, using CellListMap.
 mutable struct InteratomicDistanceFG{D,T}
     f::T
+    g::Vector{MoleculePosition{D,T}}
     dmin::T
     fmol::Vector{T} # contribution of each molecule to the function value
-    gxcar::Vector{SVector{D,T}} # Auxiliary array for gradient
+    # Auxiliary array for gradient: carries the gradient relative to the cartesian coordinates of each atom
+    gxcar::Vector{SVector{D,T}}
+    # Gradient, expressed as Vector{MoleculePosition{D,T}}, which can be reinterpted as a Vector{T}
 end
 
 # Custom copy, reset and reducer functions
 import CellListMap: copy_output, reset_output!, reducer
 function copy_output(x::InteratomicDistanceFG)
     InteratomicDistanceFG(
         x.f, 
+        copy(x.g),
         x.dmin, 
         copy(x.fmol),
         copy(x.gxcar),
     )
 end
-function reset_output!(output::InteratomicDistanceFG{N,T}) where {N,T}
+function reset_output!(output::InteratomicDistanceFG{D,T}) where {D,T}
     output.f = zero(T)
+    fill!(output.g, zero(MoleculePositions{D,T})) 
     output.dmin = typemax(T)
     fill!(output.fmol, zero(T))
-    fill!(output.gxcar, zero(SVector{N,T})) 
+    fill!(output.gxcar, zero(SVector{D,T})) 
     return output
 end
 function reducer(x::InteratomicDistanceFG, y::InteratomicDistanceFG)
     x.f += y.f
+    x.g += y.g
     x.dmin = min(x.dmin, y.dmin)
     x.fmol .+= y.fmol
     x.gxcar .+= y.gxcar
@@ -74,13 +80,11 @@ function fg!(g, system, packmol_system)
     )
     # Use the chain rule to compute the gradient relative to the rotations
     # and translations of the molecules
-    chain_rule!(packmol_system, fg)
+    chain_rule!(fg, packmol_system)
     return system.fg.f
 end
 
-function f(
-
-function pack_monoatomic_callback(spgresult, system, tol, iprint)
+function pack(spgresult, system, tol, iprint)
     if spgresult.nit % iprint == 0
         println(
             " Iteration: ", spgresult.nit,
@@ -95,39 +99,10 @@ function pack_monoatomic_callback(spgresult, system, tol, iprint)
 end
 
 """
-    pack_monoatomic!(positions::AbstractVector{<:SVector{N,T}}, unitcell, tol)
-
-Pack a monoatomic system with iniital positions `x` and distance tolerance `tol`,
-into the unitcell defined by `unitcell`, considering periodic boundary conditions.
-
-The unitcell can be a vector, in the case of orthorhombic cells, or a matrix, in the
-case of triclinic cells.
-
-The coordinates and the unitcells can be two- or three-dimensional. 
-
-# Example
-
-```julia-repl
-julia> using Packmol, StaticArrays
-
-julia> coordinates = 100 * rand(SVector{3,Float64}, 10^5);
-
-julia> unitcell = [100.0, 100.0, 100.0]; tolerance = 2.0;
-
-julia> pack_monoatomic!(coordinates, unitcell, tolerance)
-```
-
-After packing, the `coordinates` array will have updated positions for the
-atoms without overlaps. 
+    packmol()
 
 """
-function pack_monoatomic!(
-    positions::AbstractVector{<:SVector{N,T}},
-    unitcell,
-    tol;
-    parallel::Bool=true,
-    iprint=10
-) where {N,T}
+function packmol()
     #  The gradient vector will be allocated by SPGBox, as an auxiliary array
     x = copy(reinterpret(reshape, T, positions))
     auxvecs = SPGBox.VAux(x, zero(T))

src/rigid_body/chain_rule.jl

@@ -6,19 +6,24 @@ the gradient computed from the displacement of the atoms in cartesian coordinate
 
 This was implemented originally in the Fortran code of Packmol, by J. M. Martínez.
 
-Extended here for the 2D case.
+>> The function modifies the fg.g vector
 
 =#
-function chain_rule!(packmol_system::PackmolSystem{D,T}, fg) where {D,T}
-    fill!(zero(MoleculePositions{D,T}), packmol_system.gradient)
+function chain_rule!(fg, packmol_system::PackmolSystem{D,T}) where {D,T}
     imol = 0
     iat = 0
     for structure_type in packmol_system.structure_types
         for _ in 1:structure_type.number_of_molecules
             imol += 1
             ifmol = iat + 1
             ilmol = iat + structure_type.natoms
-            partial_derivatives!(imol, structure_type, packmol_system, @view(fg.gxcar[ifmol:ilmol]))
+            gmol = partial_derivatives(
+                fg.g[imol],
+                packmol_system.molecule_positon[imol].angles,
+                structure_type.reference_coordinates, 
+                @view(fg.gxcar[ifmol:ilmol]),
+            )
+            fg.g[imol] = gmol
             iat += structure_type.natoms
         end
     end
@@ -27,13 +32,21 @@ end
 
 #=
 
-3D case
+For a single molecule, computes the partial derivatives of the gradient of the
+minimum distance between atoms relative to the rotations and translations of the
+molecule, given the gradient of the minimum distance between atoms relative to the
+displacement of the atoms in cartesian coordinates.
 
 =#
-function partial_derivatives!(packmol_system, imol::Int, structure_type, gxcar)
-    g_cm = packmol_system.gradient[imol].cm # vector of gradient relative to the center of mass
-    g_rot = packmol_system.gradient[imol].angles # vector of gradient relative to the rotations
-    (sb, cb), (sg, cg), (st, ct) = sincos.(packmol_system.molecule_position[imol].angles)
+function partial_derivatives(
+    gmol::MoleculePosition{D,T}, 
+    angles::SVector{D,T},
+    reference_coordinates::Vector{SVector{D,T}},
+    gxcar::SVector{D,T},
+) where {D,T}
+    gcm = gmol.cm
+    grot = gmol.angles
+    (sb, cb), (sg, cg), (st, ct) = sincos.(angles)
     #!format off
     ∂v∂β = sum(SMatrix[
         -cb*sg*ct-sb*cg  -cb*cg*ct+sb*sg     cb*st
@@ -52,14 +65,13 @@ function partial_derivatives!(packmol_system, imol::Int, structure_type, gxcar)
     ]; dims=1)
     ∂rot = vcat(∂v∂β, ∂v∂γ, ∂v∂θ)
     #!format on
-    for i in eachindex(structure_type.reference_coordinates, gxcar)
-        x = structure_type.reference_coordinates[i]
+    for i in eachindex(reference_coordinates, gxcar)
+        x = reference_coordinates[i]
         gx = gxcar[i]
-        g_cm += gx
-        g_rot += x' * ∂rot * gx
-        packmol_system.gradient = MoleculePositions{D,T}(g_cm, g_rot)
+        gcm += gx
+        grot += x' * ∂rot * gx
     end
-    return packmol_system.gradient
+    return MoleculePosition{D,T}(gcm, grot)
 end
 
 #=