Skip to content

Commit

Permalink
Fix crystal reshaping bug
Browse files Browse the repository at this point in the history
The error was a comparison of `isapprox` for two numbers close to zero.
The new logic is simpler and more direct.
  • Loading branch information
kbarros committed Sep 28, 2023
1 parent 6494d4d commit 08e0b7c
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 16 deletions.
13 changes: 6 additions & 7 deletions src/Symmetry/Crystal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -440,7 +440,7 @@ function reshape_crystal(cryst::Crystal, new_cell_shape::Mat3)
# Lattice vectors of the new unit cell in global coordinates
new_latvecs = origin.latvecs * new_cell_shape

# Return original crystal if possible
# Return this crystal if possible
new_latvecs cryst.latvecs && return cryst

# Symmetry precision needs to be rescaled for the new unit cell. Ideally we
Expand All @@ -466,7 +466,7 @@ function reshape_crystal(cryst::Crystal, new_cell_shape::Mat3)
# Fortunately, it is well protected against mistakes via the check on atom
# count. Any mistake will yield an assertion error: "Missing atoms in
# reshaped unit cell".
nmax = sum.(eachrow(abs.(inv(B)))) .+ 1
nmax = round.(Int, sum.(eachrow(abs.(inv(B))))) .+ 1

new_positions = Vec3[]
new_types = String[]
Expand All @@ -479,11 +479,10 @@ function reshape_crystal(cryst::Crystal, new_cell_shape::Mat3)
y = B*x

# Check whether the new position y (in fractional coordinates
# associated with `new_latvecs`) is within the new unit cell. Is
# each component of y within the range [0,1)? The check
# `wrap_to_unit_cell(y) == y` accounts for finite precision ϵ by
# wrapping components of `y` to the range [-ϵ,1-ϵ).
if wrap_to_unit_cell(y; symprec=new_symprec) y
# associated with `new_latvecs`) is within the new unit cell.
# Account for finite symprec ϵ by checking the bounds [-ϵ,1-ϵ). See
# related comment in `wrap_to_unit_cell`.
if all(-new_symprec .<= y .< 1 - new_symprec)
push!(new_positions, y)
push!(new_types, cryst.types[i])
push!(new_classes, cryst.classes[i])
Expand Down
18 changes: 9 additions & 9 deletions src/Symmetry/SymmetryAnalysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,22 +8,22 @@ function all_integer(x; symprec)
return norm(x - round.(x)) < symprec
end

function is_periodic_copy(cryst::Crystal, r1::Vec3, r2::Vec3)
all_integer(r1-r2; cryst.symprec)
function is_periodic_copy(r1::Vec3, r2::Vec3; symprec)
all_integer(r1-r2; symprec)
end

function is_periodic_copy(cryst::Crystal, b1::BondPos, b2::BondPos)
function is_periodic_copy(b1::BondPos, b2::BondPos; symprec)
# Displacements between the two bonds
D1 = b2.ri - b1.ri
D2 = b2.rj - b1.rj
# Round components of D1 to nearest integers
n = round.(D1, RoundNearest)
# If both n ≈ D1 and n ≈ D2, then the bonds are equivalent by translation
return norm(n - D1) < cryst.symprec && norm(n - D2) < cryst.symprec
return norm(n - D1) < symprec && norm(n - D2) < symprec
end

function position_to_atom(cryst::Crystal, r::Vec3)
return findfirst(r′ -> is_periodic_copy(cryst, r, r′), cryst.positions)
return findfirst(r′ -> is_periodic_copy(r, r′; cryst.symprec), cryst.positions)
end

function position_to_atom_and_offset(cryst::Crystal, r::Vec3)
Expand All @@ -41,7 +41,7 @@ function symmetries_for_pointgroup_of_atom(cryst::Crystal, i::Int)
r = cryst.positions[i]
for s in cryst.symops
r′ = transform(s, r)
if is_periodic_copy(cryst, r, r′)
if is_periodic_copy(r, r′; cryst.symprec)
push!(ret, s)
end
end
Expand All @@ -56,7 +56,7 @@ function symmetries_between_atoms(cryst::Crystal, i1::Int, i2::Int)
r1 = cryst.positions[i1]
r2 = cryst.positions[i2]
for s in cryst.symops
if is_periodic_copy(cryst, r1, transform(s, r2))
if is_periodic_copy(r1, transform(s, r2); cryst.symprec)
push!(ret, s)
end
end
Expand Down Expand Up @@ -133,9 +133,9 @@ function symmetries_between_bonds(cryst::Crystal, b1::BondPos, b2::BondPos)
ret = Tuple{SymOp, Bool}[]
for s in cryst.symops
b2′ = transform(s, b2)
if is_periodic_copy(cryst, b1, b2′)
if is_periodic_copy(b1, b2′; cryst.symprec)
push!(ret, (s, true))
elseif is_periodic_copy(cryst, b1, reverse(b2′))
elseif is_periodic_copy(b1, reverse(b2′); cryst.symprec)
push!(ret, (s, false))
end
end
Expand Down

2 comments on commit 08e0b7c

@kbarros
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/92454

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.5.5 -m "<description of version>" 08e0b7cd290b01bdc1051f7d72cd92ea2b55433b
git push origin v0.5.5

Please sign in to comment.