From 193859f5619cbeca6bad9502137a800e00c63ce9 Mon Sep 17 00:00:00 2001 From: Leandro Martinez Date: Tue, 8 Oct 2024 15:40:17 -0300 Subject: [PATCH] move gradient to fg structure --- src/data_structures/PackmolSystem.jl | 1 - src/interatomic_distance_fg.jl | 51 +++++++--------------------- src/rigid_body/chain_rule.jl | 42 +++++++++++++++-------- 3 files changed, 40 insertions(+), 54 deletions(-) diff --git a/src/data_structures/PackmolSystem.jl b/src/data_structures/PackmolSystem.jl index 5d61366..8bf7060 100644 --- a/src/data_structures/PackmolSystem.jl +++ b/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) diff --git a/src/interatomic_distance_fg.jl b/src/interatomic_distance_fg.jl index 0b6fcff..6a2cbef 100644 --- a/src/interatomic_distance_fg.jl +++ b/src/interatomic_distance_fg.jl @@ -5,9 +5,12 @@ 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 @@ -15,20 +18,23 @@ 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)) diff --git a/src/rigid_body/chain_rule.jl b/src/rigid_body/chain_rule.jl index 96ba6a8..e8fdda6 100644 --- a/src/rigid_body/chain_rule.jl +++ b/src/rigid_body/chain_rule.jl @@ -6,11 +6,10 @@ 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 @@ -18,7 +17,13 @@ function chain_rule!(packmol_system::PackmolSystem{D,T}, fg) where {D,T} 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 #=