v0.3.4 (#96)

* angle calculation tweaks

* bond angle calc

* primitive temp file

* symmetry rule bug fix

* symmetry rule bug fix

Project.toml

@@ -1,7 +1,7 @@
 name = "Xtals"
 uuid = "ede5f01d-793e-4c47-9885-c447d1f18d6d"
 authors = ["SimonEnsemble <[email protected]>"]
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 Bio3DView = "99c8bb3a-9d13-5280-9740-b4880ed9c598"

docs/src/bonds.md

@@ -96,7 +96,7 @@ To get the angle (in radians) between two bonds, use `bond_angle`:
 ```jldoctest bonds
 bond_angle(xtal, 1, 5, 9)
 # output
-2.089762039489443
+2.089762039489374
 ```
 
 ## Detailed Docs

docs/src/crystal.md

@@ -24,7 +24,7 @@ xtal.bonds                          # Bonding information in the structure. By d
 xtal.symmetry                       # Symmetry information of the crystal. By default converts the symmetry to P1 symmetry.
                                     #  Use `convert_to_p1=false` argument in `Crystal` to keep original symmetry
 # output
-
+Xtals.SymmetryInfo(["x"; "y"; "z"], "P1", true)
 ```
 
 ## Fixing atom species labels
@@ -98,7 +98,18 @@ other_charged_xtal = Crystal(xtal.name, xtal.box, xtal.atoms,               # He
                              new_charges, xtal.bonds, xtal.symmetry)        #   The number of charges in the array has to be equal to the number of atoms
                                                                             #   and finally a new `Crystal` object is manually created
 # output
-
+Name: IRMOF-1.cif
+Bravais unit cell of a crystal.
+	Unit cell angles α = 90.000000 deg. β = 90.000000 deg. γ = 90.000000 deg.
+	Unit cell dimensions a = 25.832000 Å. b = 25.832000 Å, c = 25.832000 Å
+	Volume of unit cell: 17237.492730 ų
+
+	# atoms = 424
+	# charges = 424
+	chemical formula: Dict(:Zn => 4, :H => 12, :O => 13, :C => 24)
+	space Group: P1
+	symmetry Operations:
+		'x, y, z'
 ```
 
 ## Writing crystal files

src/bonds.jl

@@ -584,18 +584,21 @@ Otherwise, returns `NaN`
 """
 function bond_angle(bonds::MetaGraph, i::Int, j::Int, k::Int)::Float64
     # θ = arccos(u⃗•v⃗ / (‖u⃗‖*‖v⃗‖))
-    @assert has_prop(bonds, :has_vectors)
-    if has_edge(bonds, i, j) && has_edge(bonds, j, k)
-        # get vectors in correct orientation
-        u⃗ = get_bond_vector(bonds, j, i)
-        v⃗ = get_bond_vector(bonds, j, k)
-        norm_u⃗ = get_prop(bonds, i, j, :distance)
-        norm_v⃗ = get_prop(bonds, j, k, :distance)
-        # rounding to 12 digits avoids the case of arccos(ϕ) : ϕ > 1
-        return acos(round(dot(u⃗, v⃗) / (norm_u⃗ * norm_v⃗), digits=12))
-    else
-        return NaN
+    @assert has_prop(bonds, :has_vectors) "Vectors have not been calculated. Do so with `calculate_bond_vectors!`"
+    @assert has_edge(bonds, i, j) && has_edge(bonds, j, k) "Invalid vertex selection. Cannot calculate ∠($i,$j,$k)"
+    # get vectors in correct orientation
+    u⃗ = get_bond_vector(bonds, j, i)
+    v⃗ = get_bond_vector(bonds, j, k)
+    norm_u⃗ = get_prop(bonds, i, j, :distance)
+    norm_v⃗ = get_prop(bonds, j, k, :distance)
+    ϕ = dot(u⃗, v⃗) / (norm_u⃗ * norm_v⃗)
+    # correct for numerical precision domain error |ϕ| > 1
+    if isapprox(ϕ, 1)
+        ϕ = 1.0
+    elseif isapprox(ϕ, -1)
+        ϕ = -1.0
     end
+    return acos(ϕ)
 end
 
 bond_angle(xtal::Crystal, i::Int, j::Int, k::Int) = bond_angle(xtal.bonds, i, j, k)

src/crystal.jl

@@ -809,6 +809,9 @@ function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=t
 
     @printf(cif_file, "_symmetry_Int_Tables_number 1\n\n")
     @printf(cif_file, "loop_\n_symmetry_equiv_pos_as_xyz\n")
+    if size(crystal.symmetry.operations, 2) == 0
+        @printf(cif_file, "'x,y,z'\n")
+    end
     for i in 1:size(crystal.symmetry.operations, 2)
         @printf(cif_file, "'%s,%s,%s'\n", crystal.symmetry.operations[:, i]...)
     end
@@ -1106,7 +1109,7 @@ function primitive_cell(xtal::Crystal)
     else
         pymatgen = rc[:pymatgen]
     end
-    tempfile = rc[:paths][:crystals] * "/temp_$(uuid1()).cif"
+    tempfile = joinpath(pwd(), ".temp_$(uuid1()).cif")
     # copy out xtal and convert it to primitive cell
     write_cif(xtal, tempfile)
     pymatgen.CifParser(tempfile).get_structures()[1].get_primitive_structure().to(filename=tempfile)

test/bonds.jl

@@ -52,7 +52,7 @@ end
     # check bond angle
     @test isapprox(bond_angle(xtal, 1, 2, 3), deg2rad(90))
     # check invalid bond angle
-    @test isnan(bond_angle(xtal, 2, 3, 1))
+    @test_throws AssertionError bond_angle(xtal, 2, 3, 1)
     # check bond distance calculations
     @test isapprox(bond_distance(xtal, 1, 2), norm(get_prop(xtal.bonds, 1, 2, :vector)))
 

test/crystal.jl

@@ -23,9 +23,9 @@ end
     @test Xtals.vtk_filename(xtal) == "ATIBOU01_clean.vtk"
 
     # primitive cells via pymatgen
-    xtal = Crystal("IRMOF-1.cif")
+    xtal = Crystal("str_m1_o3_o20_pcu_sym.22.cif")
     prim = primitive_cell(xtal)
-    @test prim.atoms.n == 106
+    @test prim.atoms.n == 996
     @test isnothing(assert_P1_symmetry(prim))
     pymatgen = rc[:pymatgen]
     rc[:pymatgen] = nothing