Skip to content

Commit

Permalink
bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@MatthiasSachs
MatthiasSachs committed Oct 15, 2024
1 parent a8d78fe commit bda49bd
Show file tree
Hide file tree
Showing 3 changed files with 94 additions and 28 deletions.
16 changes: 14 additions & 2 deletions src/frictionfit/fluxmodels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,15 +36,27 @@ function _Gamma(B::AbstractArray{T,5}, cc::AbstractArray{T,2}, ::Type{<:ACMatrix
return Γ
end

# function _Gamma(Bt::AbstractArray{T,4}, cc::AbstractArray{T,2}, ::Type{<:PWCMatrixModel}) where {T}
# @tullio Σ[d,i,j,r] := Bt[d,i,j,k] * cc[k,r]
# @tullio Γ[d1,d2,i,j] := - Σ[d1,i,j,r] * Σ[d2,j,i,r]
# @tullio Γd[d1,d2,i,i] := Σ[d1,i,j,r] * Σ[d2,i,j,r]
# return Γ + Γd
# end
# function _Gamma(Bt::AbstractArray{T,5}, cc::AbstractArray{T,2}, ::Type{<:PWCMatrixModel}) where {T}
# @tullio Σ[d1,d2,i,j,r] := Bt[d1,d2,i,j,k] * cc[k,r]
# @tullio Γ[d1,d2,i,j] := - Σ[d1,d,i,j,r] * Σ[d2,d,j,i,r]
# @tullio Γd[d1,d2,i,i] := Σ[d1,d,i,j,r] * Σ[d2,d,i,j,r]
# return Γ + Γd
# end
function _Gamma(Bt::AbstractArray{T,4}, cc::AbstractArray{T,2}, ::Type{<:PWCMatrixModel}) where {T}
@tullio Σ[d,i,j,r] := Bt[d,i,j,k] * cc[k,r]
@tullio Γ[d1,d2,i,j] := - Σ[d1,i,j,r] * Σ[d2,j,i,r]
@tullio Γ[d1,d2,i,j] := Σ[d1,i,j,r] * Σ[d2,j,i,r]
@tullio Γd[d1,d2,i,i] := Σ[d1,i,j,r] * Σ[d2,i,j,r]
return Γ + Γd
end
function _Gamma(Bt::AbstractArray{T,5}, cc::AbstractArray{T,2}, ::Type{<:PWCMatrixModel}) where {T}
@tullio Σ[d1,d2,i,j,r] := Bt[d1,d2,i,j,k] * cc[k,r]
@tullio Γ[d1,d2,i,j] := - Σ[d1,d,i,j,r] * Σ[d2,d,j,i,r]
@tullio Γ[d1,d2,i,j] := Σ[d1,d,i,j,r] * Σ[d2,d,j,i,r]
@tullio Γd[d1,d2,i,i] := Σ[d1,d,i,j,r] * Σ[d2,d,i,j,r]
return Γ + Γd
end
Expand Down
64 changes: 47 additions & 17 deletions src/frictionmodels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,32 +95,68 @@ function Gamma(M::PWCMatrixModel, at::Atoms; kvargs...)
return sum(_square(Σ,M) for Σ in Σ_vec)
end

# function _square(Σ::SparseMatrixCSC{Tv,Ti}, ::PWCMatrixModel) where {Tv, Ti}
# # if !iseven(length(Σ.nzval))
# # @warn "Self-interactions of particles are present and included friction matrix."
# # end
# #@assert iseven(length(Σ.nzval))
# nvals = 2*length(Σ.nzval) #+ length(Σ.m)
# Is, Js, Vs = findnz(Σ)
# I, J, V = Ti[],Ti[],SMatrix{3, 3,eltype(Tv), 9}[]
# #Array{Ti}(undef,nvals), Array{Ti}(undef,nvals), Array{SMatrix{3, 3,eltype(Tv), 9}}(undef,nvals)
# sizehint!(I, nvals)
# sizehint!(J, nvals)
# sizehint!(V, nvals)
# #k = 1
# for (i,j,σij) in zip(Is, Js, Vs)
# if i <= j
# # if i == j
# # @warn "Self-interactions of particles are present and included friction matrix."
# # end
# σji = Σ[j,i]
# push!(I, i)
# push!(J,j)
# push!(V, -σij * σji')

# push!(I, j)
# push!(J, i)
# push!(V, -σji * σij')

# push!(I, i)
# push!(J, i)
# push!(V, σij* σij')

# push!(I, j)
# push!(J, j)
# push!(V, σji* σji')
# #I[k], J[k], V[k] = i,j,-σij * σji'
# #I[k+1], J[k+1], V[k+1] = j,i, -σji * σij'
# #I[k+2], J[k+2], V[k+2] = i,i, σij* σij'
# #I[k+3], J[k+3], V[k+3] = j,j, σji* σji'
# #k+=4
# end
# end
# A = sparse(I, J, V, Σ.m, Σ.n)
# return A
# end

function _square(Σ::SparseMatrixCSC{Tv,Ti}, ::PWCMatrixModel) where {Tv, Ti}
# if !iseven(length(Σ.nzval))
# @warn "Self-interactions of particles are present and included friction matrix."
# end
#@assert iseven(length(Σ.nzval))
nvals = 2*length(Σ.nzval) #+ length(Σ.m)
Is, Js, Vs = findnz(Σ)
I, J, V = Ti[],Ti[],SMatrix{3, 3,eltype(Tv), 9}[]
#Array{Ti}(undef,nvals), Array{Ti}(undef,nvals), Array{SMatrix{3, 3,eltype(Tv), 9}}(undef,nvals)
sizehint!(I, nvals)
sizehint!(J, nvals)
sizehint!(V, nvals)
#k = 1
for (i,j,σij) in zip(Is, Js, Vs)
if i <= j
# if i == j
# @warn "Self-interactions of particles are present and included friction matrix."
# end
σji = Σ[j,i]
push!(I, i)
push!(J,j)
push!(V, -σij * σji')
push!(V, σij * σji')

push!(I, j)
push!(J, i)
push!(V, -σji * σij')
push!(V, σji * σij')

push!(I, i)
push!(J, i)
Expand All @@ -129,19 +165,13 @@ function _square(Σ::SparseMatrixCSC{Tv,Ti}, ::PWCMatrixModel) where {Tv, Ti}
push!(I, j)
push!(J, j)
push!(V, σji* σji')
#I[k], J[k], V[k] = i,j,-σij * σji'
#I[k+1], J[k+1], V[k+1] = j,i, -σji * σij'
#I[k+2], J[k+2], V[k+2] = i,i, σij* σij'
#I[k+3], J[k+3], V[k+3] = j,j, σji* σji'
#k+=4
end
end
A = sparse(I, J, V, Σ.m, Σ.n)
return A
end



# using Tullio
# function Gamma(M::MatrixModel{Covariant}, at::Atoms; kvargs...)
# Σ_vec = Sigma(M, at; kvargs...)
Expand Down
42 changes: 33 additions & 9 deletions src/matrixmodels/matrixmodels.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,14 @@ abstract type SpeciesCoupling end
struct SpeciesCoupled <: SpeciesCoupling end
struct SpeciesUnCoupled <: SpeciesCoupling end

function ACE.write_dict(sc::SC) where {SC<:SpeciesCoupling}
return Dict("__id__" => string("ACEds_SpeciesCoupling"), "sc"=>typeof(sc))
end
function ACE.read_dict(::Val{:ACEds_SpeciesCoupling}, D::Dict)
sc = getfield(ACEds.MatrixModels, Symbol(D["sc"]))
return sc()
end

abstract type NoiseCoupling end

struct PairCoupling <: NoiseCoupling end
Expand All @@ -81,20 +89,34 @@ function ACE.read_dict(::Val{:ACEds_NoiseCoupling}, D::Dict)
return coupling()
end

# _mreduce(z1::Symbol,z2::Symbol, sc::SpeciesCoupled) = chemical_symbol.(_mreduce(AtomicNumber(z1),AtomicNumber(z2), sc))
# _mreduce(z1::Symbol,z2::Symbol, ::Type{SpeciesCoupled}) = chemical_symbol.(_mreduce(AtomicNumber(z1),AtomicNumber(z2), SpeciesCoupled))
_mreduce(z1,z2, ::SpeciesUnCoupled) = (z1,z2)
_mreduce(z1,z2, ::SpeciesCoupled) = _msort(z1,z2)
_mreduce(z1,z2, ::Type{SpeciesUnCoupled}) = (z1,z2)
_mreduce(z1,z2, ::Type{SpeciesCoupled}) = _msort(z1,z2)

#TODO: needs to be corrected because the function will falsely return SpeciesCoupled() if only one species
function _species_symmetry(mkeys)
if all([(z2,z1)==_msort(z1,z2) for (z1,z2) in mkeys])
return SpeciesCoupled()
elseif all([(z2,z1) in mkeys for (z1,z2) in mkeys])
return SpeciesUnCoupled()
else
@error "The species coupling is inconsistent."
end
function _assert_consistency(mkeys, ::SpeciesUnCoupled)
return @assert all([((z2,z1) in mkeys && (z1,z2) in mkeys) for (z1,z2) in mkeys])
end

function _assert_consistency(mkeys, ::SpeciesCoupled)
return @assert all([ begin (z1s,z2s) = _msort(z1,z2);
((z1s,z2s) in mkeys && ((z1s==z2s) || !((z2s,z1s) in mkeys)))
end for (z1,z2) in mkeys])
end

#TODO: needs to be corrected because the function will falsely return SpeciesCoupled() if only one species
# function _species_symmetry(mkeys) # this function is not working!!
# if all([(z2,z1)==_msort(z1,z2) for (z1,z2) in mkeys])
# return SpeciesCoupled()
# elseif all([(z2,z1) in mkeys for (z1,z2) in mkeys])
# return SpeciesUnCoupled()
# else
# @error "The species coupling is inconsistent."
# end
# end

function _assert_offsite_keys(offsite_dict, ::SpeciesCoupled)
return @assert all([(z2,z1)==_msort(z1,z2) for (z1,z2) in keys(offsite_dict)])
end
Expand All @@ -112,6 +134,8 @@ _T(::ACE.LinearACEModel{TB, SVector{N,T}, TEV}) where {TB,N,T,TEV} = T


_msort(z1,z2) = (z1<=z2 ? (z1,z2) : (z2,z1))
# TODO: it may be better to base sorting on Atomic numbers instead of chemical_symbols
_msort(z1::AtomicNumber,z2::AtomicNumber) = map(AtomicNumber,_msort(chemical_symbol(z1),chemical_symbol(z2)))

NamedCollection = Union{AbstractDict,NamedTuple}

Expand Down

0 comments on commit bda49bd

Please sign in to comment.