From f99aa603e8b0ef158a14476173c7805c2b1edcb8 Mon Sep 17 00:00:00 2001 From: Adrian Henle Date: Wed, 12 Oct 2022 11:25:51 -0700 Subject: [PATCH] formatting (#178) * formatting * JuliaFormatter Action Bot (#179) Co-authored-by: eahenle * JuliaFormatter Action Bot (#180) Co-authored-by: eahenle Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: eahenle --- .github/workflow_extras/write_load_combos.jl | 6 +- README.md | 4 +- docs/make.jl | 36 +- docs/src/bonds.md | 44 +- docs/src/box.md | 12 +- docs/src/crystal.md | 42 +- docs/src/distance.md | 8 +- docs/src/globals.md | 20 +- docs/src/index.md | 12 +- docs/src/matter.md | 58 +- docs/src/visualization.md | 2 +- quick_setup.jl | 11 +- src/Xtals.jl | 88 ++- src/bonds.jl | 299 +++++--- src/box.jl | 146 ++-- src/crystal.jl | 686 +++++++++++------ src/distance.jl | 98 +-- src/matter.jl | 65 +- src/misc.jl | 120 +-- src/repfactors.jl | 27 +- test/aqua.jl | 7 +- test/assert_p1_symmetry.jl | 2 +- test/bond_vectors.jl | 92 +-- test/bonds.jl | 64 +- test/box.jl | 92 +-- test/crystal.jl | 747 ++++++++++--------- test/distance.jl | 80 +- test/matter.jl | 125 ++-- test/misc.jl | 19 +- test/paths.jl | 2 +- test/runtests.jl | 6 +- 31 files changed, 1802 insertions(+), 1218 deletions(-) diff --git a/.github/workflow_extras/write_load_combos.jl b/.github/workflow_extras/write_load_combos.jl index b24167b..9854afd 100644 --- a/.github/workflow_extras/write_load_combos.jl +++ b/.github/workflow_extras/write_load_combos.jl @@ -1,7 +1,7 @@ using Combinatorics open("load_combos.sh", "w") do f - for (a,b,c,d) in permutations(ARGS) + for (a, b, c, d) in permutations(ARGS) write(f, "julia --project -e 'using $a; using $b; using $c; using $d' && \\\n") end - write(f, "echo -e \\\\nDone!\n") -end \ No newline at end of file + return write(f, "echo -e \\\\nDone!\n") +end diff --git a/README.md b/README.md index a649413..c1ebaa2 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ ![Xtals.jl](logo.jpg) -| **Documentation** | **Build Status** | **Test Coverage** | -|:---:|:---:|:---:| +| **Documentation** | **Build Status** | **Test Coverage** | +|:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------:|:-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------:|:--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------:| | [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://SimonEnsemble.github.io/Xtals.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://SimonEnsemble.github.io/Xtals.jl/dev) | [![CI](https://github.com/SimonEnsemble/Xtals.jl/actions/workflows/CI_Build.yml/badge.svg)](https://github.com/SimonEnsemble/Xtals.jl/actions/workflows/CI_Build.yml) [![Docs](https://github.com/SimonEnsemble/Xtals.jl/actions/workflows/doc_deployment.yml/badge.svg)](https://github.com/SimonEnsemble/Xtals.jl/actions/workflows/doc_deployment.yml) [![weekly](https://github.com/SimonEnsemble/Xtals.jl/actions/workflows/weekly.yml/badge.svg)](https://github.com/SimonEnsemble/Xtals.jl/actions/workflows/weekly.yml) | [![codecov](https://codecov.io/gh/SimonEnsemble/Xtals.jl/branch/master/graph/badge.svg?token=QM6XZ3KAW1)](https://codecov.io/gh/SimonEnsemble/Xtals.jl) [![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) | A pure-[Julia](https://julialang.org/) package for importing, manipulating, and working with crystal structures, such as metal-organic frameworks (MOFs). diff --git a/docs/make.jl b/docs/make.jl index 017ab0d..96a2639 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,23 +1,23 @@ using Documenter using Xtals, XtalsPyTools -makedocs( - root = joinpath(dirname(pathof(Xtals)), "..", "docs"), - modules = [Xtals, XtalsPyTools], - sitename = "Xtals.jl", - clean = true, - pages = [ - "Xtals" => "index.md", - "globals" => "globals.md", - "matter" => "matter.md", - "boxes" => "box.md", - "crystals" => "crystal.md", - "computing distances" => "distance.md", - "bonding" => "bonds.md", - "visualization" => "visualization.md", - "XtalsPyTools" => "pytools.md" - ], - format = Documenter.HTML(assets = ["assets/flux.css"]) +makedocs(; + root=joinpath(dirname(pathof(Xtals)), "..", "docs"), + modules=[Xtals, XtalsPyTools], + sitename="Xtals.jl", + clean=true, + pages=[ + "Xtals" => "index.md", + "globals" => "globals.md", + "matter" => "matter.md", + "boxes" => "box.md", + "crystals" => "crystal.md", + "computing distances" => "distance.md", + "bonding" => "bonds.md", + "visualization" => "visualization.md", + "XtalsPyTools" => "pytools.md" + ], + format=Documenter.HTML(; assets=["assets/flux.css"]) ) -deploydocs(repo = "github.com/SimonEnsemble/Xtals.jl.git") +deploydocs(; repo="github.com/SimonEnsemble/Xtals.jl.git") diff --git a/docs/src/bonds.md b/docs/src/bonds.md index 4276ad6..2f92215 100644 --- a/docs/src/bonds.md +++ b/docs/src/bonds.md @@ -11,34 +11,38 @@ The nodes of the graph correspond to the [`Crystal`](@ref)'s [`Atoms`](@ref), an ## Bonding Rules -`Xtals` uses an array of [`BondingRule`](@ref) structs stored in `rc` for deciding if two [`Atoms`](@ref) are an appropriate distance to be chemically bonded. -The default rules are based on the [Cordero covalent radii](https://doi.org/10.1039/B801115J), modified based on the work of [Thomas Manz](https://doi.org/10.1039/c9ra07327b). +`Xtals` uses an array of [`BondingRule`](@ref) structs stored in `rc` for deciding if two [`Atoms`](@ref) are an appropriate distance to be chemically bonded. +The default rules are based on the [Cordero covalent radii](https://doi.org/10.1039/B801115J), modified based on the work of [Thomas Manz](https://doi.org/10.1039/c9ra07327b). Each [`BondingRule`](@ref) is composed of two chemical species symbols and a floating point value, the maximum distance for inferring a bond between the indicated species. ```jldoctest; output=false BondingRule(:C, :C, 1.77) + # output + BondingRule(:C, :C, 1.77) ``` -The global bonding rules may be augmented with [`add_bonding_rules`](@ref) or written to/read from disk with [`write_bonding_rules`](@ref) and [`read_bonding_rules`](@ref). -The default rules are determined from `rc[:covalent_radii]` at module load, but *are not recalculated upon changes to the covalent radii.* +The global bonding rules may be augmented with [`add_bonding_rules`](@ref) or written to/read from disk with [`write_bonding_rules`](@ref) and [`read_bonding_rules`](@ref). +The default rules are determined from `rc[:covalent_radii]` at module load, but *are not recalculated upon changes to the covalent radii.* If `rc[:covalent_radii]` is altered and new bonding rules should be calculated, the user must do `rc[:bonding_rules] = bondingrules()`. ## Adding Bonds to Crystals -By default, bonding information is not added to a [`Crystal`](@ref). -Bonds may be inferred at the time of loading structure data by using the `infer_bonds` keyword argument. +By default, bonding information is not added to a [`Crystal`](@ref). +Bonds may be inferred at the time of loading structure data by using the `infer_bonds` keyword argument. See [`Crystal`](@ref) for more details. ```jldoctest bonds -xtal = Crystal("SBMOF-1.cif", infer_bonds=true, periodic_boundaries=true) +xtal = Crystal("SBMOF-1.cif"; infer_bonds=true, periodic_boundaries=true) xtal.bonds + # output + {120, 144} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0) ``` -The first number is the number of nodes, and the second is the number of edges. +The first number is the number of nodes, and the second is the number of edges. The `:weight` attribute is not used, and can be ignored. [`remove_bonds!`](@ref) clears the bonding information from a [`Crystal`](@ref): @@ -46,7 +50,9 @@ The `:weight` attribute is not used, and can be ignored. ```jldoctest bonds remove_bonds!(xtal) xtal.bonds + # output + {120, 0} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0) ``` @@ -55,46 +61,54 @@ Use [`infer_bonds!`](@ref) to infer plausible bonds using the global bonding rul ```jldoctest bonds infer_bonds!(xtal, true) xtal.bonds + # output + {120, 144} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0) ``` ## Bonds from Input File -Some chemical information formats, like `.cif` and `.mol`, can store not only the cartesian coordinates of atoms, but also the graph of bonds between atoms in molecules and crystals. -The `read_bonds_from_file` keyword argument for [`Crystal`](@ref) enables loading these bonds when reading the data. +Some chemical information formats, like `.cif` and `.mol`, can store not only the cartesian coordinates of atoms, but also the graph of bonds between atoms in molecules and crystals. +The `read_bonds_from_file` keyword argument for [`Crystal`](@ref) enables loading these bonds when reading the data. [`read_mol`](@ref) also returns bond information. ## Bonds for Atoms -In the case that atomic coordinates are loaded from XYZ format, there will be no unit cell information. +In the case that atomic coordinates are loaded from XYZ format, there will be no unit cell information. To infer bonds between atoms in this case, use the [`infer_bonds`](@ref) function: ```jldoctest bonds atoms = Cart(xtal.atoms, xtal.box) # get atoms in Cartesian coords bonds = infer_bonds(atoms) # infer bonding graph + # output + {120, 110} undirected Int64 metagraph with Float64 weights defined by :weight (default weight 1.0) ``` ## Bond Distances, Vectors, and Angles -Bonds may be labeled with several additional pieces of information. +Bonds may be labeled with several additional pieces of information. The first is the center-to-center distance between the bonded atoms, accessible via [`bond_distance`](@ref): ```jldoctest bonds bond_distance(xtal, 1, 5) + # output + 1.5233240952030063 ``` -The bond distance is automatically added by [`infer_bonds!`](@ref). +The bond distance is automatically added by [`infer_bonds!`](@ref). Applying [`calculate_bond_vectors!`](@ref) (or passing `calculate_vectors=true` to [`infer_bonds!`](@ref)) labels each edge in the bonding graph with a vector, accessible via [`get_bond_vector`](@ref): ```jldoctest bonds calculate_bond_vectors!(xtal) get_bond_vector(xtal, 1, 5) + # output + 3-element Vector{Float64}: -1.2460713575602618 -0.26441824999999985 @@ -105,7 +119,9 @@ While the bond graph itself is undirected, the vectors returned by [`get_bond_ve ```jldoctest bonds get_bond_vector(xtal, 5, 1) + # output + 3-element Vector{Float64}: 1.2460713575602618 0.26441824999999985 @@ -117,7 +133,9 @@ To get the angle (in radians) between two bonds, use [`bond_angle`](@ref): ```jldoctest bonds bond_angle(xtal, 1, 5, 9) + # output + 2.089762039489374 ``` diff --git a/docs/src/box.md b/docs/src/box.md index 69648fd..02a44be 100644 --- a/docs/src/box.md +++ b/docs/src/box.md @@ -18,9 +18,9 @@ For example, given the unit cell of Co-MOF-74, we can define its [`Box`](@ref): a = 26.13173 # Å b = 26.13173 c = 6.722028 -α = π/2 # radians -β = π/2 -γ = 2*π/3 +α = π / 2 # radians +β = π / 2 +γ = 2 * π / 3 box = Box(a, b, c, α, β, γ) ``` @@ -40,7 +40,6 @@ To quickly get a simple unit-cubic [`Box`](@ref), use the [`unit_cube`](@ref) fu #└ Volume of unit cell: 1.000000 ų ``` - ## Transforming Coordinates Conversions are provided for switching between [`Frac`](@ref)tional and [`Cart`](@ref)esian [`Coords`](@ref) using the [`Box`](@ref) (works for [`Atoms`](@ref) and [`Charges`](@ref), too). @@ -53,14 +52,13 @@ Cart(xtal.atoms.coords, xtal.box) # 1.231811631 0.32198514120000005 … 6.2082409932000004 2.2119953472]) ``` - ## Replicating a Box For simulations in larger volumes than a single crystallograhic unit cell, the [`Box`](@ref) may be replicated along each or any of the three crystallographic axes. See [`replicate`](@ref). ```julia -replicated_box = replicate(box, (2,2,2)) +replicated_box = replicate(box, (2, 2, 2)) ``` ## Exporting a Box @@ -69,7 +67,9 @@ For visualization of the unit cell boundaries, the [`Box`](@ref) may be written ```jldoctest; setup=:(box = Box([26.1317 -13.0659 0; 0 22.6307 0; 0 0 6.72203])), output=false write_vtk(box, "box.vtk"; verbose=true) + # output + See box.vtk ``` diff --git a/docs/src/crystal.md b/docs/src/crystal.md index 2976d84..fa6747f 100644 --- a/docs/src/crystal.md +++ b/docs/src/crystal.md @@ -10,8 +10,8 @@ end ## Reading in a Crystal Structure File -Currently, the crystal structure file reader accepts `.cif` and `.cssr` file formats. -`Xtals.jl` looks for the crystal structure files in `rc[:paths][:crystals]` which is by default `./data/crystals` (relative to `pwd()` at module loading). +Currently, the crystal structure file reader accepts `.cif` and `.cssr` file formats. +`Xtals.jl` looks for the crystal structure files in `rc[:paths][:crystals]` which is by default `./data/crystals` (relative to `pwd()` at module loading). By typing `rc[:paths][:crystals] = "my_crystal_dir"`, `Xtals.jl` now looks for the crystal structure file in `my_crystal_dir`. The files can be read as: @@ -22,37 +22,43 @@ xtal.box # The unit cell information xtal.atoms # The atom coordinates (in fractional space) and the atom identities xtal.charges # The charge magnitude and coordinates (in fractional space) xtal.bonds # Bonding information in the structure. By default this is an empty graph, - # but use `read_bonds_from_file=true` argument in `Crystal` to read from crystal structure file +# but use `read_bonds_from_file=true` argument in `Crystal` to read from crystal structure file 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 +# Use `convert_to_p1=false` argument in `Crystal` to keep original symmetry + # output + Xtals.SymmetryInfo(["x"; "y"; "z";;], "P1", true) ``` ## Fixing Atom Species Labels -Often, the atoms species are appended by numbers. +Often, the atoms species are appended by numbers. This messes with the internal workings of `Xtals.jl`. To circumvent this problem, the function [`strip_numbers_from_atom_labels!`](@ref) removes the appended numbers. It is important to use this function prior to GCMC or Henry coefficient calculations and bond inference operations. ```jldoctest crystal -xtal = Crystal("IRMOF-1.cif", species_col=["_atom_site_label"]) +xtal = Crystal("IRMOF-1.cif"; species_col=["_atom_site_label"]) xtal.atoms.species[1] + # output + :Zn1 ``` ```jldoctest crystal strip_numbers_from_atom_labels!(xtal) xtal.atoms.species[1] + # output + :Zn ``` ## Converting the Coordinates to Cartesian Space -The coordinates of the [`Crystal`](@ref)'s [`Atoms`](@ref) are stored in [`Frac`](@ref)tional coordinates. +The coordinates of the [`Crystal`](@ref)'s [`Atoms`](@ref) are stored in [`Frac`](@ref)tional coordinates. If one needs to analyze the [`Cart`](@ref)esian coordinates of the [`Crystal`](@ref), that can be done by using the unit cell ([`Box`](@ref)) information. ```julia @@ -66,9 +72,11 @@ For many simulations, one needs to replicate the unit cell multiple times to cre This is done with [`replicate`](@ref): ```jldoctest crystal -super_xtal = replicate(xtal, (2,2,2)) # Replicates the original unit cell once in each dimension +super_xtal = replicate(xtal, (2, 2, 2)) # Replicates the original unit cell once in each dimension xtal.atoms.n, super_xtal.atoms.n + # output + (424, 3392) ``` @@ -87,12 +95,20 @@ If the structure file does not contain partial charges, we provide methods to as ```jldoctest crystal; output=false species_to_charge = Dict(:Zn => 2.0, :C => 0.0, :H => 0.0, :O => -0.61538) # This method assigns a static charge to atom species charged_xtal = assign_charges(xtal, species_to_charge, 1e-3) # This function creates a new charged `Crystal` object. - # The function checks for charge neutrality with a tolerance of 1e-3 +# The function checks for charge neutrality with a tolerance of 1e-3 new_charges = Charges(zeros(xtal.atoms.n), xtal.atoms.coords) -other_charged_xtal = Crystal(xtal.name, xtal.box, xtal.atoms, # Here we create a new `Charges` object using an array of new charges. - 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 +other_charged_xtal = Crystal( + xtal.name, + xtal.box, + xtal.atoms, # Here we create a new `Charges` object using an array of new charges. + 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. @@ -117,11 +133,11 @@ We provide methods to write both `.xyz` and `.cif` files. ```jldoctest crystal; output=false write_cif(xtal, "my_new_cif_file.cif") # Stored in the current directory write_xyz(xtal, "my_new_xyz_file.xyz") # stored in the current directory + # output ``` - # Detailed Docs ```@docs diff --git a/docs/src/distance.md b/docs/src/distance.md index e048897..e9b02b4 100644 --- a/docs/src/distance.md +++ b/docs/src/distance.md @@ -6,13 +6,15 @@ end # Distances -The distance between two [`Atoms`](@ref) in a [`Crystal`](@ref) is central to many operations within `Xtals.jl`. +The distance between two [`Atoms`](@ref) in a [`Crystal`](@ref) is central to many operations within `Xtals.jl`. The [`distance`](@ref) function calculates the [`Cart`](@ref)esian displacement between the [`Coords`](@ref) ([`Cart`](@ref) or [`Frac`](@ref)) of two points, `i` and `j`, within a given [`Box`](@ref), in units of Å. ```jldoctest distance xtal = Crystal("IRMOF-1.cif") distance(xtal.atoms.coords, xtal.box, 1, 2, false) + # output + 18.538930020700434 ``` @@ -20,7 +22,9 @@ The `apply_pbc` argument allows for calculation of distances across the periodic ```jldoctest distance distance(xtal.atoms.coords, xtal.box, 1, 2, true) + # output + 15.096469110488975 ``` @@ -28,7 +32,9 @@ distance(xtal.atoms.coords, xtal.box, 1, 2, true) ```jldoctest distance distance(xtal.atoms, xtal.box, 3, 5, true) + # output + 16.90555095103936 ``` diff --git a/docs/src/globals.md b/docs/src/globals.md index 5fd3c43..a2aa62e 100644 --- a/docs/src/globals.md +++ b/docs/src/globals.md @@ -1,23 +1,23 @@ # Global Variables -In `Xtals.jl`, global variables are stored in a dictionary called `rc`. +In `Xtals.jl`, global variables are stored in a dictionary called `rc`. The entries in `rc` are used for holding various important pieces of information: -- `rc[:atomic_masses]` : a dictionary which maps atomic species to their masses in AMU -- `rc[:cpk_colors]` : a dictionary which maps atomic species to their CPK colors -- `rc[:bonding_rules]` : an array which lists the maximum bonding distances between each pair of atom types in Å -- `rc[:covalent_radii]` : a dictionary which maps atomic species to their covalent radii in Å -- `rc[:paths]` : a dictionary of paths to directories to search for input data + - `rc[:atomic_masses]` : a dictionary which maps atomic species to their masses in AMU + - `rc[:cpk_colors]` : a dictionary which maps atomic species to their CPK colors + - `rc[:bonding_rules]` : an array which lists the maximum bonding distances between each pair of atom types in Å + - `rc[:covalent_radii]` : a dictionary which maps atomic species to their covalent radii in Å + - `rc[:paths]` : a dictionary of paths to directories to search for input data ## Paths The `rc[:paths]` dictionary holds the following: -- `rc[:paths][:crystals]` : the path to find crystallographic files (`.cif`, `.cssr`) to read in -- `rc[:paths][:data]` : not directly used by `Xtals.jl`, but important for other packages which use `Xtals`, e.g. [PorousMaterials.jl](https://simonensemble.github.io/PorousMaterials.jl/dev/) + - `rc[:paths][:crystals]` : the path to find crystallographic files (`.cif`, `.cssr`) to read in + - `rc[:paths][:data]` : not directly used by `Xtals.jl`, but important for other packages which use `Xtals`, e.g. [PorousMaterials.jl](https://simonensemble.github.io/PorousMaterials.jl/dev/) -When `Xtals` is first loaded, the paths are set based on the present working directory. -`Xtals` expects to find a folder named `data` in the working directory, and assumes that crystallographic data will be present in `data/crystals`. +When `Xtals` is first loaded, the paths are set based on the present working directory. +`Xtals` expects to find a folder named `data` in the working directory, and assumes that crystallographic data will be present in `data/crystals`. To change the location of the crystal inputs, either set `rc[:paths][:crystals]` directly, or use [`set_paths`](@ref): ```julia diff --git a/docs/src/index.md b/docs/src/index.md index 3e89e32..cc38b94 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -6,15 +6,15 @@ A pure-[Julia](https://julialang.org/) package for representation of porous crys ## Installation -1. Download and install the [Julia programming language](https://julialang.org/), v1.3 or higher. + 1. Download and install the [Julia programming language](https://julialang.org/), v1.3 or higher. -2. Open the package manager in the Julia REPL (using `]`) and enter: + 2. Open the package manager in the Julia REPL (using `]`) and enter: ``` pkg> add Xtals ``` -3. In Julia, load all functions in `Xtals.jl` into the namespace: + 3. In Julia, load all functions in `Xtals.jl` into the namespace: ``` julia> using Xtals # that's it @@ -22,9 +22,9 @@ julia> using Xtals # that's it ## Dependencies -Some (non-core) features of `Xtals.jl` require Python libraries. -If these are missing, you may see warnings when the module loads. -For users unfamiliar with configuring Python, the [quick-setup script](https://raw.githubusercontent.com/SimonEnsemble/Xtals.jl/master/quick_setup.jl) will install and build `Xtals.jl` along with a Python environment and the dependencies. +Some (non-core) features of `Xtals.jl` require Python libraries. +If these are missing, you may see warnings when the module loads. +For users unfamiliar with configuring Python, the [quick-setup script](https://raw.githubusercontent.com/SimonEnsemble/Xtals.jl/master/quick_setup.jl) will install and build `Xtals.jl` along with a Python environment and the dependencies. Run it like so: ``` diff --git a/docs/src/matter.md b/docs/src/matter.md index be55adc..5864bc5 100644 --- a/docs/src/matter.md +++ b/docs/src/matter.md @@ -6,12 +6,12 @@ end # Matter and Coordinates -[`Atoms`](@ref) and [`Charges`](@ref) are the building blocks of [`Crystal`](@ref)s and molecules in `Xtals.jl`. +[`Atoms`](@ref) and [`Charges`](@ref) are the building blocks of [`Crystal`](@ref)s and molecules in `Xtals.jl`. Each have coordinates in both [`Cart`](@ref)esian and [`Frac`](@ref)tional space (associated with unit cell information, i.e., a [`Box`](@ref)). ## Coordinates -We store coordinates as an abstract [`Coords`](@ref) type that has two subtypes: `Cart` and `Frac` for Cartesian and Fractional, respectively. +We store coordinates as an abstract [`Coords`](@ref) type that has two subtypes: `Cart` and `Frac` for Cartesian and Fractional, respectively. See the [Wikipedia](https://en.wikipedia.org/wiki/Fractional_coordinates) page on fractional coordinates, which are defined in the context of a periodic system, e.g. within a crystal. Construct coordinates of `n` particles by passing a `n` by `3` array: @@ -20,7 +20,9 @@ Construct coordinates of `n` particles by passing a `n` by `3` array: # construct cartesian coordinates of a particle coord = Cart([1.0, 2.0, 5.0]) coord.x + # output + 3×1 Matrix{Float64}: 1.0 2.0 @@ -31,7 +33,9 @@ coord.x # construct fractional coordinates of a particle coord = Frac([0.1, 0.2, 0.5]) coord.xf + # output + 3×1 Matrix{Float64}: 0.1 0.2 @@ -43,11 +47,13 @@ The coordinates of multiple particles are stored column-wise: ```jldoctest matter; output=false # five particles at uniform random coordinates coords = Cart([ - 0.0 1.0 0.0 0.0 1.0; - 0.0 0.0 1.0 0.0 1.0; - 0.0 0.0 0.0 1.0 1.0 + 0.0 1.0 0.0 0.0 1.0 + 0.0 0.0 1.0 0.0 1.0 + 0.0 0.0 0.0 1.0 1.0 ]) + # output + Cart([0.0 1.0 … 0.0 1.0; 0.0 0.0 … 0.0 1.0; 0.0 0.0 … 1.0 1.0]) ``` @@ -59,7 +65,9 @@ coords[2:3] # (slicing by index) coords of particles 2 and 3 coords[[1, 2, 5]] # (slicing by index) coords of particles 1, 2, and 5 coords[rand(Bool, 5)] # (boolean slicing) coords, selected at random length(coords) # number of particles, (5) + # output + 5 ``` @@ -69,7 +77,9 @@ length(coords) # number of particles, (5) ```jldoctest matter coords.x = rand(3, 5) + # output + ERROR: setfield!: immutable struct of type Cart cannot be changed ``` @@ -86,7 +96,9 @@ Fractional coordinates can be wrapped to be inside the unit cell box: coords = Frac([1.2, -0.3, 0.9]) wrap!(coords) coords.xf + # output + 3×1 Matrix{Float64}: 0.19999999999999996 0.7 @@ -97,10 +109,12 @@ We can translate coordinates by a vector `dx`: ```jldoctest dx = Cart([1.0, 2.0, 3.0]) -coords = Cart([1.0, 0.0, 0.0]) +coords = Cart([1.0, 0.0, 0.0]) translate_by!(coords, dx) coords.x + # output + 3×1 Matrix{Float64}: 2.0 2.0 @@ -115,7 +129,9 @@ box = unit_cube() coords = Cart([1.0, 0.0, 0.0]) translate_by!(coords, dx, box) coords.x + # output + 3×1 Matrix{Float64}: 1.1 0.20000000000000004 @@ -128,16 +144,19 @@ An atom is specified by its coordinates and atomic species. We can construct a s ```jldoctest matter; output=false species = [:O, :H, :H] # atomic species are represnted with Symbols -coords = Cart([0.0 0.757 -0.757; # coordinates of each - 0.0 0.586 0.586; - 0.0 0.0 0.0 ] - ) +coords = Cart([ + 0.0 0.757 -0.757 # coordinates of each + 0.0 0.586 0.586 + 0.0 0.0 0.0 +]) atoms = Atoms(species, coords) # 3 atoms comprising water atoms.n # number of atoms, 3 atoms.coords # coordinates; atoms.coords.x gives the array of coords atoms.species # array of species atoms::Atoms{Cart} # successful type assertion, as opposed to atoms::Atoms{Frac} + # output + Atoms{Cart}(3, [:O, :H, :H], Cart([0.0 0.757 -0.757; 0.0 0.586 0.586; 0.0 0.0 0.0])) ``` @@ -148,7 +167,9 @@ We can slice [`Atoms`](@ref), such as: ```jldoctest matter; output=false atoms[1] # 1st atom atoms[2:3] # 2nd and 3rd atom + # output + Atoms{Cart}(2, [:H, :H], Cart([0.757 -0.757; 0.586 0.586; 0.0 0.0])) ``` @@ -157,7 +178,9 @@ And combine them: ```jldoctest matter atoms_combined = atoms[1] + atoms[2:3] # combine atoms 1, 2, and 3 isapprox(atoms, atoms_combined) + # output + true ``` @@ -167,16 +190,19 @@ Point [`Charges`](@ref) work analogously to [`Atoms`](@ref), except instead of ` ```jldoctest matter; output=false q = [-1.0, 0.5, 0.5] # values of point charges, units: electrons -coords = Cart([0.0 0.757 -0.757; # coordinates of the point charges - 0.0 0.586 0.586; - 0.0 0.0 0.0 ] - ) +coords = Cart([ + 0.0 0.757 -0.757 # coordinates of the point charges + 0.0 0.586 0.586 + 0.0 0.0 0.0 +]) charges = Charges(q, coords) # 3 point charges charges.n # number of charges, 3 charges.coords # retreive coords charges.q # retreive q charges::Charges{Cart} # successful type assertion, as opposed to charges::Charges{Frac} + # output + Charges{Cart}(3, [-1.0, 0.5, 0.5], Cart([0.0 0.757 -0.757; 0.0 0.586 0.586; 0.0 0.0 0.0])) ``` @@ -184,13 +210,17 @@ We can determine if the set of point charges comprise a charge-neutral system by ```jldoctest matter net_charge(charges) + # output + 0.0 ``` ```jldoctest matter neutral(charges) + # output + true ``` diff --git a/docs/src/visualization.md b/docs/src/visualization.md index 2b94038..e6f9b5f 100644 --- a/docs/src/visualization.md +++ b/docs/src/visualization.md @@ -2,7 +2,7 @@ Crystals can be visualized using [`Bio3DView`](https://github.com/jgreener64/Bio3DView.jl) and either [`Blink`](https://juliagizmos.github.io/Blink.jl/stable/), [`Jupyter`](https://julialang.github.io/IJulia.jl/stable/), or [`Pluto`](https://github.com/fonsp/Pluto.jl). -To include a visualization within a `Pluto` or `Jupyter` notebook, simply call [`view_crystal`](@ref). +To include a visualization within a `Pluto` or `Jupyter` notebook, simply call [`view_crystal`](@ref). To launch a visualization window directly from a script or the REPL, use `Blink`, like so: ```julia diff --git a/quick_setup.jl b/quick_setup.jl index aa18745..3192bce 100644 --- a/quick_setup.jl +++ b/quick_setup.jl @@ -1,22 +1,19 @@ import Pkg - # tries to load a Python dependency; on failure, adds the dependency via Conda -function check_add_dep(pkg; channel="") +check_add_dep(pkg; channel="") = try @info "Checking dependency: $pkg" pyimport(pkg) catch @info "Installing $pkg..." - Conda.add(pkg, channel=channel) + Conda.add(pkg; channel=channel) pyimport(pkg) @info "$pkg install verified." end -end - @info "Setting up Python environment..." -ENV["PYTHON"]="" +ENV["PYTHON"] = "" Pkg.add("PyCall") Pkg.build("PyCall") using PyCall @@ -63,6 +60,6 @@ end # check the deps, add if missing check_add_dep("scipy") -check_add_dep("pymatgen", channel="conda-forge") +check_add_dep("pymatgen"; channel="conda-forge") @info "Setup complete!" diff --git a/src/Xtals.jl b/src/Xtals.jl index 4ab7a57..5aa3fb4 100644 --- a/src/Xtals.jl +++ b/src/Xtals.jl @@ -1,11 +1,21 @@ module Xtals -using AtomsBase, Bio3DView, CSV, DataFrames, FIGlet, Graphs, LinearAlgebra, MetaGraphs, Printf, StaticArrays, Unitful +using AtomsBase, + Bio3DView, + CSV, + DataFrames, + FIGlet, + Graphs, + LinearAlgebra, + MetaGraphs, + Printf, + StaticArrays, + Unitful using PrecompileSignatures: @precompile_signatures # global variable dictionary -global rc = Dict{Symbol,Any}() -rc[:paths] = Dict{Symbol,String}() +global rc = Dict{Symbol, Any}() +rc[:paths] = Dict{Symbol, String}() include("matter.jl") include("box.jl") @@ -24,42 +34,88 @@ rc[:covalent_radii] = DEFAULT_COVALENT_RADII DEFAULT_BONDING_RULES = bondingrules() rc[:bonding_rules] = DEFAULT_BONDING_RULES - function __init__() # create path entries in global dictionary rc[:paths][:crystals] = "" rc[:paths][:data] = "" # set paths to data and crystals relative to pwd() at import - set_paths(joinpath(pwd(), "data"); no_warn=true) + return set_paths(joinpath(pwd(), "data"); no_warn=true) end - export # Xtals.jl rc, # matter.jl - Coords, Frac, Cart, Atoms, Charges, wrap!, neutral, net_charge, translate_by!, + Coords, + Frac, + Cart, + Atoms, + Charges, + wrap!, + neutral, + net_charge, + translate_by!, # box.jl - Box, replicate, unit_cube, write_vtk, inside, + Box, + replicate, + unit_cube, + write_vtk, + inside, # distance.jl - nearest_image!, distance, overlap, remove_duplicates, pairwise_distances, + nearest_image!, + distance, + overlap, + remove_duplicates, + pairwise_distances, # misc.jl - read_xyz, write_xyz, read_mol, write_mol2, assert_P1_symmetry, set_paths, view_crystal, + read_xyz, + write_xyz, + read_mol, + write_mol2, + assert_P1_symmetry, + set_paths, + view_crystal, # crystal.jl - Crystal, strip_numbers_from_atom_labels!, assign_charges, empirical_formula, molecular_weight, - crystal_density, write_cif, has_charges, apply_symmetry_operations, write_cssr, rename_xtal, + Crystal, + strip_numbers_from_atom_labels!, + assign_charges, + empirical_formula, + molecular_weight, + crystal_density, + write_cif, + has_charges, + apply_symmetry_operations, + write_cssr, + rename_xtal, # AtomsBase things also from crystal - position, velocity, bounding_box, boundary_conditions, chemical_formula, + position, + velocity, + bounding_box, + boundary_conditions, + chemical_formula, # bonds.jl - infer_bonds!, write_bond_information, BondingRule, bond_sanity_check, remove_bonds!, read_bonding_rules, - write_bonding_rules, add_bonding_rules, drop_cross_pb_bonds!, bondingrules, bond_angle, - get_bond_vector, calculate_bond_vectors!, clear_vectors!, bond_distance, infer_bonds + infer_bonds!, + write_bond_information, + BondingRule, + bond_sanity_check, + remove_bonds!, + read_bonding_rules, + write_bonding_rules, + add_bonding_rules, + drop_cross_pb_bonds!, + bondingrules, + bond_angle, + get_bond_vector, + calculate_bond_vectors!, + clear_vectors!, + bond_distance, + infer_bonds @precompile_signatures(Xtals) diff --git a/src/bonds.jl b/src/bonds.jl index 91697cc..832a428 100644 --- a/src/bonds.jl +++ b/src/bonds.jl @@ -6,9 +6,10 @@ A rule for determining if two atoms within a crystal are bonded. # Attributes -- `species_i::Symbol`: One of the atoms types for this bond rule -- `species_j::Symbol`: The other atom type for this bond rule -- `max_dist`: The maximum distance between the atoms for bonding to occur + + - `species_i::Symbol`: One of the atoms types for this bond rule + - `species_j::Symbol`: The other atom type for this bond rule + - `max_dist`: The maximum distance between the atoms for bonding to occur """ struct BondingRule species_i::Symbol @@ -16,9 +17,9 @@ struct BondingRule max_dist::Float64 end -BondingRule(species_i::String, species_j::String, max_dist::String) = BondingRule( - Symbol.([species_i, species_j])..., parse(Float64, max_dist)) - +function BondingRule(species_i::String, species_j::String, max_dist::String) + return BondingRule(Symbol.([species_i, species_j])..., parse(Float64, max_dist)) +end """ bond_rules = bondingrules() @@ -28,13 +29,17 @@ Returns a set of bonding rules based on the given Cordero parameters and toleran # Arguments -- `cov_rad::Union{Dict{Symbol, Float64}, Nothing}`: Covalent radii. See [`covalent_radii()`](@ref) -- `pad::Float`: Amount to pad covalent radii for bonding interactions. + - `cov_rad::Union{Dict{Symbol, Float64}, Nothing}`: Covalent radii. See [`covalent_radii()`](@ref) + - `pad::Float`: Amount to pad covalent radii for bonding interactions. # Returns -- `bondingrules::Array{BondingRule, 1}`: Bonding rules from the specified covalent radii.` + + - `bondingrules::Array{BondingRule, 1}`: Bonding rules from the specified covalent radii.` """ -function bondingrules(;covalent_radii::Dict{Symbol,Float64}=rc[:covalent_radii], pad::Float64=0.)::Array{BondingRule} +function bondingrules(; + covalent_radii::Dict{Symbol, Float64}=rc[:covalent_radii], + pad::Float64=0.0 +)::Array{BondingRule} bond_rules = BondingRule[] # loop over parameterized atoms for (i, atom_i) in enumerate(keys(covalent_radii)) @@ -51,29 +56,38 @@ function bondingrules(;covalent_radii::Dict{Symbol,Float64}=rc[:covalent_radii], return bond_rules end - """ write_bonding_rules("file.csv") + Writes bonding rules to a CSV file that can be loaded with [`read_bonding_rules`](@ref) + # Arguments -- `filename::String` : The name of the output file -- `bonding_rules::Array{BondingRule}` : (Optional) The rules to write to file. If not specified, the global rules are written. + + - `filename::String` : The name of the output file + - `bonding_rules::Array{BondingRule}` : (Optional) The rules to write to file. If not specified, the global rules are written. """ -function write_bonding_rules(filename::String, bonding_rules::Array{BondingRule}=rc[:bonding_rules]) +function write_bonding_rules( + filename::String, + bonding_rules::Array{BondingRule}=rc[:bonding_rules] +) open(filename, "w") do f - for r ∈ bonding_rules + for r in bonding_rules @printf(f, "%s,%s,%f\n", r.species_i, r.species_j, r.max_dist) end end end - """ read_bonding_rules("file.csv") + Reads a CSV file of bonding rules and returns a BondingRule array. + # Arguments -- `filename::String` : name of file in data directory containing bonding rules + + - `filename::String` : name of file in data directory containing bonding rules + # Returns + `rules::Array{BondingRule}` : the bonding rules read from file """ function read_bonding_rules(filename::String)::Array{BondingRule} @@ -86,20 +100,21 @@ function read_bonding_rules(filename::String)::Array{BondingRule} return rules end - """ add_bonding_rules(bonding_rules) + Adds `bonding_rules` to the beginning of the global bonding rules list + # Arguments -- `bonding_rules::Array{BondingRule}` : the array of bonding rules to add + + - `bonding_rules::Array{BondingRule}` : the array of bonding rules to add """ function add_bonding_rules(bonding_rules::Array{BondingRule}) br = rc[:bonding_rules] - rc[:bonding_rules] = vcat(bonding_rules, br) + return rc[:bonding_rules] = vcat(bonding_rules, br) end add_bonding_rules(bonding_rule::BondingRule) = add_bonding_rules([bonding_rule]) - # for pretty-printing the bonding rules function show(io::IO, bonding_rules::Array{BondingRule}) for r in bonding_rules @@ -113,18 +128,25 @@ end Checks to see if atoms `i` and `j` in `crystal` are bonded according to the `bonding_rules`. # Arguments -- `atoms::Atoms{Cart}`: The atoms that may be bonded -- `i::Int`: Index of the first atom -- `j::Int`: Index of the second atom -- `bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will + + - `atoms::Atoms{Cart}`: The atoms that may be bonded + - `i::Int`: Index of the first atom + - `j::Int`: Index of the second atom + - `bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will be used to fill the bonding information. They are applied in the order that they appear. # Returns -- `bonded_flag::Bool`: true iff atoms `i` and `j` are bonded according to `bonding_rules` -- `r::Float64`: the distance between atomic centers. + + - `bonded_flag::Bool`: true iff atoms `i` and `j` are bonded according to `bonding_rules` + - `r::Float64`: the distance between atomic centers. """ -function is_bonded(atoms::Atoms{Cart}, i::Int64, j::Int64, bonding_rules::Array{BondingRule, 1}) +function is_bonded( + atoms::Atoms{Cart}, + i::Int64, + j::Int64, + bonding_rules::Array{BondingRule, 1} +) species_i = atoms.species[i] species_j = atoms.species[j] # compute vector between the two atoms. @@ -140,18 +162,26 @@ end """ bonded = loop_over_bonding_rules(bonding_rules, species_i, species_j, r) """ -function loop_over_bonding_rules(bonding_rules::Vector{BondingRule}, species_i::Symbol, species_j::Symbol, r::Float64) +function loop_over_bonding_rules( + bonding_rules::Vector{BondingRule}, + species_i::Symbol, + species_j::Symbol, + r::Float64 +) # loop over possible bonding rules for br in bonding_rules # determine if the atom species correspond to the species in `bonding_rules` species_match = false if br.species_i == :* && br.species_j == :* species_match = true - elseif br.species_i == :* && (species_i == br.species_j || species_j == br.species_j) + elseif br.species_i == :* && + (species_i == br.species_j || species_j == br.species_j) species_match = true - elseif br.species_j == :* && (species_i == br.species_i || species_j == br.species_j) + elseif br.species_j == :* && + (species_i == br.species_i || species_j == br.species_j) species_match = true - elseif (species_i == br.species_i && species_j == br.species_j) || (species_j == br.species_i && species_i == br.species_j) + elseif (species_i == br.species_i && species_j == br.species_j) || + (species_j == br.species_i && species_i == br.species_j) species_match = true end if species_match @@ -169,7 +199,6 @@ function loop_over_bonding_rules(bonding_rules::Vector{BondingRule}, species_i:: return false end - """ are_atoms_bonded = is_bonded(crystal, i, j, bonding_rules=[BondingRule(:H, :*, 1.2), BondingRule(:*, :*, 1.9)], include_bonds_across_periodic_boundaries=true) @@ -177,22 +206,29 @@ end Checks to see if atoms `i` and `j` in `crystal` are bonded according to the `bonding_rules`. # Arguments -- `crystal::Crystal`: The crystal that bonds will be added to -- `i::Int`: Index of the first atom -- `j::Int`: Index of the second atom -- `bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will + + - `crystal::Crystal`: The crystal that bonds will be added to + - `i::Int`: Index of the first atom + - `j::Int`: Index of the second atom + - `bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will be used to fill the bonding information. They are applied in the order that they appear. -- `include_bonds_across_periodic_boundaries::Bool`: Whether to check across the + - `include_bonds_across_periodic_boundaries::Bool`: Whether to check across the periodic boundary when calculating bonds # Returns -- `bonded_flag::Bool`: true iff atoms `i` and `j` are bonded according to `bonding_rules` -- `r::Float64`: the periodic distance between atomic centers. For non-bonded pairs, `r=NaN`. -- `cross_boundary_flag::Bool`: true iff the bond is across a periodic boundary. -""" -function is_bonded(crystal::Crystal, i::Int64, j::Int64, bonding_rules::Array{BondingRule, 1}; - include_bonds_across_periodic_boundaries::Bool=true) + + - `bonded_flag::Bool`: true iff atoms `i` and `j` are bonded according to `bonding_rules` + - `r::Float64`: the periodic distance between atomic centers. For non-bonded pairs, `r=NaN`. + - `cross_boundary_flag::Bool`: true iff the bond is across a periodic boundary. +""" +function is_bonded( + crystal::Crystal, + i::Int64, + j::Int64, + bonding_rules::Array{BondingRule, 1}; + include_bonds_across_periodic_boundaries::Bool=true +) species_i = crystal.atoms.species[i] species_j = crystal.atoms.species[j] # compute vector between the two atoms. @@ -218,7 +254,6 @@ function is_bonded(crystal::Crystal, i::Int64, j::Int64, bonding_rules::Array{Bo return bonded, r, cross_boundary_flag end - """ remove_bonds!(crystal) @@ -226,28 +261,31 @@ Remove all bonds from a crystal structure, `crystal::Crystal`. """ function remove_bonds!(crystal::Crystal) while ne(crystal.bonds) > 0 - rem_edge!(crystal.bonds, collect(edges(crystal.bonds))[1].src, - collect(edges(crystal.bonds))[1].dst) + rem_edge!( + crystal.bonds, + collect(edges(crystal.bonds))[1].src, + collect(edges(crystal.bonds))[1].dst + ) end - clear_vectors!(crystal) + return clear_vectors!(crystal) end - """ bonds = infer_bonds(atoms; bonding_rules=rc[:bonding_rules]) Returns a `MetaGraph` encoding the chemical bonds between atoms. # Arguments -- `atoms::Atoms{Cart}`: the atoms to bond -- `bonding_rules::Vector{BondingRule}`: the bonding rules to use + + - `atoms::Atoms{Cart}`: the atoms to bond + - `bonding_rules::Vector{BondingRule}`: the bonding rules to use """ function infer_bonds(atoms::Atoms{Cart}; bonding_rules=rc[:bonding_rules])::MetaGraph bonds = MetaGraph(atoms.n) - for i in 1:atoms.n + for i in 1:(atoms.n) # loop over every unique pair of atoms - for j in i+1:atoms.n + for j in (i + 1):(atoms.n) bonding_flag, r = is_bonded(atoms, i, j, bonding_rules) if bonding_flag add_edge!(bonds, i, j) @@ -260,7 +298,6 @@ function infer_bonds(atoms::Atoms{Cart}; bonding_rules=rc[:bonding_rules])::Meta return bonds end - """ infer_bonds!(crystal, include_bonds_across_periodic_boundaries; bonding_rules=rc[:bonding_rules]) @@ -271,21 +308,35 @@ Populate the bonds in the crystal object based on the bonding rules. If a pair d The bonding rules are hierarchical, i.e. the first bonding rule takes precedence over the latter ones. # Arguments -- `crystal::Crystal`: The crystal that bonds will be added to. -- `include_bonds_across_periodic_boundaries::Bool`: Whether to check across the periodic boundary when calculating bonds. -- `bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will be used to fill the bonding information. They are applied in the order that they appear. `rc[:bonding_rules]` will be used if none provided. -- `calculate_vectors::Bool`: Optional. Set `true` to annotate all edges in the `bonds` graph with vector information. -- `sanity_check::Bool`: Optional. Set `false` to skip the sanity check after inferring bonds. -""" -function infer_bonds!(crystal::Crystal, include_bonds_across_periodic_boundaries::Bool; - bonding_rules::Array{BondingRule, 1}=rc[:bonding_rules], calculate_vectors::Bool=false, sanity_check::Bool=true) - @assert ne(crystal.bonds) == 0 @sprintf("The crystal %s already has bonds. Remove them with the `remove_bonds!` function before inferring new ones.", crystal.name) + + - `crystal::Crystal`: The crystal that bonds will be added to. + - `include_bonds_across_periodic_boundaries::Bool`: Whether to check across the periodic boundary when calculating bonds. + - `bonding_rules::Array{BondingRule, 1}`: The array of bonding rules that will be used to fill the bonding information. They are applied in the order that they appear. `rc[:bonding_rules]` will be used if none provided. + - `calculate_vectors::Bool`: Optional. Set `true` to annotate all edges in the `bonds` graph with vector information. + - `sanity_check::Bool`: Optional. Set `false` to skip the sanity check after inferring bonds. +""" +function infer_bonds!( + crystal::Crystal, + include_bonds_across_periodic_boundaries::Bool; + bonding_rules::Array{BondingRule, 1}=rc[:bonding_rules], + calculate_vectors::Bool=false, + sanity_check::Bool=true +) + @assert ne(crystal.bonds) == 0 @sprintf( + "The crystal %s already has bonds. Remove them with the `remove_bonds!` function before inferring new ones.", + crystal.name + ) # loop over every atom - for i in 1:crystal.atoms.n + for i in 1:(crystal.atoms.n) # loop over every unique pair of atoms - for j in i+1:crystal.atoms.n - bonding_flag, r, cross_boundary_flag = is_bonded(crystal, i, j, bonding_rules; - include_bonds_across_periodic_boundaries=include_bonds_across_periodic_boundaries) + for j in (i + 1):(crystal.atoms.n) + bonding_flag, r, cross_boundary_flag = is_bonded( + crystal, + i, + j, + bonding_rules; + include_bonds_across_periodic_boundaries=include_bonds_across_periodic_boundaries + ) if bonding_flag add_edge!(crystal.bonds, i, j) set_prop!(crystal.bonds, i, j, :distance, r) @@ -301,14 +352,14 @@ function infer_bonds!(crystal::Crystal, include_bonds_across_periodic_boundaries end end - """ sane_bonds = bond_sanity_check(crystal) Run sanity checks on `crystal.bonds`. -* is the bond graph fully connected? i.e. does every vertex (=atom) in the bond graph have at least one edge? -* each hydrogen can have only one bond -* each carbon can have a maximum of four bonds + + - is the bond graph fully connected? i.e. does every vertex (=atom) in the bond graph have at least one edge? + - each hydrogen can have only one bond + - each carbon can have a maximum of four bonds if sanity checks fail, refer to [`write_bond_information`](@ref) to write a .vtk to visualize the bonds. @@ -317,7 +368,7 @@ Return `true` if sanity checks pass, `false` otherwise. """ function bond_sanity_check(crystal::Crystal)::Bool pass_flag = true - for a = 1:crystal.atoms.n + for a in 1:(crystal.atoms.n) ns = neighbors(crystal.bonds, a) # is the graph fully connected? if length(ns) == 0 @@ -338,7 +389,6 @@ function bond_sanity_check(crystal::Crystal)::Bool return pass_flag end - """ write_bond_information(crystal) write_bond_information(crystal, filename) @@ -349,26 +399,35 @@ end Writes the bond information from a crystal to the selected filename. # Arguments -- `crystal::Crystal`: The crystal to have its bonds written to a vtk file -- `filename::String`: The filename the bond information will be saved to. If left out, will default to crystal name. -- `center_at_origin::Bool`: (optional) center the coordinates at the origin of the crystal -- `bond_filter::Pair{Symbol, Function}`: (optional) a key-value pair of an edge attribute and a predicate function. Bonds with attributes that cause the predicate to return false are excluded from writing. -- `verbose::Bool`: (optional) if true, prints output file name to console. -""" -function write_bond_information(crystal::Crystal, filename::String; - verbose::Bool=false, center_at_origin::Bool=false, bond_filter::Pair{Symbol, F}=(:NOTHING => x -> ())) where F <: Function + + - `crystal::Crystal`: The crystal to have its bonds written to a vtk file + - `filename::String`: The filename the bond information will be saved to. If left out, will default to crystal name. + - `center_at_origin::Bool`: (optional) center the coordinates at the origin of the crystal + - `bond_filter::Pair{Symbol, Function}`: (optional) a key-value pair of an edge attribute and a predicate function. Bonds with attributes that cause the predicate to return false are excluded from writing. + - `verbose::Bool`: (optional) if true, prints output file name to console. +""" +function write_bond_information( + crystal::Crystal, + filename::String; + verbose::Bool=false, + center_at_origin::Bool=false, + bond_filter::Pair{Symbol, F}=(:NOTHING => x -> ()) +) where {F <: Function} if ne(crystal.bonds) == 0 - @warn("Crystal %s has no bonds present. To get bonding information for this - crystal run `infer_bonds!` with an array of bonding rules\n", crystal.name) + @warn( + "Crystal %s has no bonds present. To get bonding information for this + crystal run `infer_bonds!` with an array of bonding rules\n", + crystal.name + ) end - if ! occursin(".vtk", filename) + if !occursin(".vtk", filename) filename *= ".vtk" end # filter bonds idx_keep_bonds = trues(ne(crystal.bonds)) attr, pred = bond_filter if attr != :NOTHING - for (b, bond) ∈ enumerate(edges(crystal.bonds)) + for (b, bond) in enumerate(edges(crystal.bonds)) prop = get_prop(crystal.bonds, bond, attr) if !pred(prop) idx_keep_bonds[b] = false @@ -377,15 +436,29 @@ function write_bond_information(crystal::Crystal, filename::String; end # write output open(filename, "w") do vtk_file - @printf(vtk_file, "# vtk DataFile Version 2.0\n%s bond information\nASCII\n - DATASET POLYDATA\nPOINTS %d double\n", crystal.name, nv(crystal.bonds)) - for i = 1:crystal.atoms.n + @printf( + vtk_file, + "# vtk DataFile Version 2.0\n%s bond information\nASCII\n +DATASET POLYDATA\nPOINTS %d double\n", + crystal.name, + nv(crystal.bonds) + ) + for i in 1:(crystal.atoms.n) if center_at_origin - @printf(vtk_file, "%0.5f\t%0.5f\t%0.5f\n", (crystal.box.f_to_c * - (crystal.atoms.coords.xf[:, i] - [0.5, 0.5, 0.5]))...) + @printf( + vtk_file, + "%0.5f\t%0.5f\t%0.5f\n", + ( + crystal.box.f_to_c * + (crystal.atoms.coords.xf[:, i] - [0.5, 0.5, 0.5]) + )... + ) else - @printf(vtk_file, "%0.5f\t%0.5f\t%0.5f\n", (crystal.box.f_to_c * - crystal.atoms.coords.xf[:, i])...) + @printf( + vtk_file, + "%0.5f\t%0.5f\t%0.5f\n", + (crystal.box.f_to_c * crystal.atoms.coords.xf[:, i])... + ) end end @printf(vtk_file, "\nLINES %d %d\n", sum(idx_keep_bonds), 3 * sum(idx_keep_bonds)) @@ -396,12 +469,17 @@ function write_bond_information(crystal::Crystal, filename::String; end end if verbose - @printf("Saved bond information for crystal %s to %s.\n", crystal.name, joinpath(pwd(), filename)) + @printf( + "Saved bond information for crystal %s to %s.\n", + crystal.name, + joinpath(pwd(), filename) + ) end end -write_bond_information(crystal::Crystal; kwargs...) = write_bond_information(crystal, crystal.name * "_bonds.vtk"; kwargs...) - +function write_bond_information(crystal::Crystal; kwargs...) + return write_bond_information(crystal, crystal.name * "_bonds.vtk"; kwargs...) +end """ Loop through xtal and calculate any missing distances @@ -411,26 +489,29 @@ function calc_missing_bond_distances!(xtal::Crystal) if !(:distance ∈ keys(props(xtal.bonds, bond))) i = src(bond) j = dst(bond) - set_prop!(xtal.bonds, i, j, :distance, distance(xtal.atoms, xtal.box, i, j, true)) + set_prop!( + xtal.bonds, + i, + j, + :distance, + distance(xtal.atoms, xtal.box, i, j, true) + ) end end end - """ Gets rid of the bonds across unit cell boundaries """ -function drop_cross_pb_bonds!(bonds::MetaGraph) +drop_cross_pb_bonds!(bonds::MetaGraph) = for bond in collect(edges(bonds)) if get_prop(bonds, bond, :cross_boundary) rem_edge!(bonds, src(bond), dst(bond)) end end -end drop_cross_pb_bonds!(xtal::Crystal) = drop_cross_pb_bonds!(xtal.bonds) - """ `clear_vectors!(xtal)` `clear_vectors!(xtal.bonds)` @@ -438,18 +519,18 @@ drop_cross_pb_bonds!(xtal::Crystal) = drop_cross_pb_bonds!(xtal.bonds) Removes edge vector attributes from crystal bonding graph. """ function clear_vectors!(bonds::MetaGraph) - for edge ∈ edges(bonds) + for edge in edges(bonds) rem_prop!(bonds, edge, :vector) rem_prop!(bonds, edge, :direction) end - set_prop!(bonds, :has_vectors, false) + return set_prop!(bonds, :has_vectors, false) end clear_vectors!(xtal::Crystal) = clear_vectors!(xtal.bonds) - """ `calculate_bond_vectors!(xtal)` + Adds a property to the edges of a crystal's bonding graph giving the vector between source and destination nodes in Cartesian space. """ function calculate_bond_vectors!(xtal::Crystal) @@ -462,13 +543,15 @@ function calculate_bond_vectors!(xtal::Crystal) # check that vectors have not already been calculated if has_prop(xtal.bonds, :has_vectors) && get_prop(xtal.bonds, :has_vectors) # loop will skip all bonds if vectors already set. fatal error. - error("Crystal $(xtal.name) already has vectors. Clear them with `clear_vectors!` first.") + error( + "Crystal $(xtal.name) already has vectors. Clear them with `clear_vectors!` first." + ) return end # get bonds as directed graph (matters for vectors!) bonds = DiGraph(adjacency_matrix(xtal.bonds)) # loop over bonding edges - for edge ∈ edges(bonds) + for edge in edges(bonds) # get node IDs i, j = src(edge), dst(edge) # only need to label an edge once in undirected graph @@ -488,12 +571,12 @@ function calculate_bond_vectors!(xtal::Crystal) # xtal.bonds is undirected, so need to specify direction of vector set_prop!(xtal.bonds, i, j, :direction, (i, j)) end - set_prop!(xtal.bonds, :has_vectors, true) + return set_prop!(xtal.bonds, :has_vectors, true) end - """ `get_bond_vector(bonds, i, j)` + Returns the vector between connected nodes `i` and `j` of the `bonds` graph. """ function get_bond_vector(bonds::MetaGraph, i::Int, j::Int) @@ -506,7 +589,6 @@ end get_bond_vector(xtal::Crystal, i::Int, j::Int) = get_bond_vector(xtal.bonds, i, j) - """ `bond_angle(xtal.bonds, 2, 1, 3)` `bond_angle(xtal, 8, 121, 42)` @@ -537,7 +619,6 @@ end bond_angle(xtal::Crystal, i::Int, j::Int, k::Int) = bond_angle(xtal.bonds, i, j, k) - """ `bond_distance(xtal, i, j)` `bond_distance(xtal.bonds, i, j)` diff --git a/src/box.jl b/src/box.jl index 2eaf87a..acc25fa 100644 --- a/src/box.jl +++ b/src/box.jl @@ -54,13 +54,20 @@ function reciprocal_lattice(f_to_c::Array{Float64, 2}) end # Automatically calculates Ω, f_to_c, and c_to_f for `Box` data structure based on axes lenghts and angles. -function Box(a::Float64, b::Float64, c::Float64, - α::Float64, β::Float64, γ::Float64) +function Box(a::Float64, b::Float64, c::Float64, α::Float64, β::Float64, γ::Float64) # unit cell volume (A³) - Ω = a * b * c * sqrt(1 - cos(α) ^ 2 - cos(β) ^ 2 - cos(γ) ^ 2 + 2 * cos(α) * cos(β) * cos(γ)) + Ω = a * b * c * sqrt(1 - cos(α)^2 - cos(β)^2 - cos(γ)^2 + 2 * cos(α) * cos(β) * cos(γ)) # matrices to map fractional coords <--> Cartesian coords - f_to_c = [[a, 0, 0] [b * cos(γ), b * sin(γ), 0] [c * cos(β), c * (cos(α) - cos(β) * cos(γ)) / sin(γ), Ω / (a * b * sin(γ))]] - c_to_f = [[1/a, 0, 0] [-cos(γ) / (a * sin(γ)), 1 / (b * sin(γ)), 0] [b * c * (cos(α) * cos(γ) - cos(β)) / (Ω * sin(γ)), a * c * (cos(β) * cos(γ) - cos(α)) / (Ω * sin(γ)), a * b * sin(γ) / Ω]] + f_to_c = [[a, 0, 0] [b * cos(γ), b * sin(γ), 0] [ + c * cos(β), + c * (cos(α) - cos(β) * cos(γ)) / sin(γ), + Ω / (a * b * sin(γ)) + ]] + c_to_f = [[1 / a, 0, 0] [-cos(γ) / (a * sin(γ)), 1 / (b * sin(γ)), 0] [ + b * c * (cos(α) * cos(γ) - cos(β)) / (Ω * sin(γ)), + a * c * (cos(β) * cos(γ) - cos(α)) / (Ω * sin(γ)), + a * b * sin(γ) / Ω + ]] # the columns of f_to_c are the unit cell axes r = reciprocal_lattice(f_to_c) @@ -101,8 +108,8 @@ end This function generates a unit cube, each side is 1.0 Angstrom long, and all the corners are right angles. """ -unit_cube() = Box(1.0, 1.0, 1.0, π/2, π/2, π/2) -Box(a::Float64, b::Float64, c::Float64) = Box(a, b, c, π/2, π/2, π/2) # right angle box +unit_cube() = Box(1.0, 1.0, 1.0, π / 2, π / 2, π / 2) +Box(a::Float64, b::Float64, c::Float64) = Box(a, b, c, π / 2, π / 2, π / 2) # right angle box """ f_coords = Frac(c_coords, box) @@ -111,7 +118,6 @@ Box(a::Float64, b::Float64, c::Float64) = Box(a, b, c, π/2, π/2, π/2) # right convert Cartesian coordinates `c_coords::Cart` to fractional coordinates `f_coords::Frac`, using the `box::Box`. this works on `Atoms` and `Charges` too. - """ Frac(coords::Cart, box::Box) = Frac(box.c_to_f * coords.x) Frac(atoms::Atoms{Cart}, box::Box) = Atoms(atoms.species, Frac(atoms.coords, box)) @@ -122,7 +128,6 @@ Frac(charges::Charges{Cart}, box::Box) = Charges(charges.q, Frac(charges.coords, atoms_c = Cart(atoms_f, box) charges_c = Cart(charges_f, box) - convert fractional coordinates `f_coords::Frac` to Cartesian coordinates `c_coords::Cart`, using the `box::Box`. this works on `Atoms` and `Charges` too. """ @@ -139,15 +144,23 @@ Note `replicate(original_box, repfactors=(1, 1, 1))` returns same `Box`. The new fractional coordinates as described by `f_to_c` and `c_to_f` still ∈ [0, 1]. # Arguments -- `original_box::Box`: The box that you want to replicate -- `repfactors::Tuple{Int, Int, Int}`: The factor you want to replicate the box by + + - `original_box::Box`: The box that you want to replicate + - `repfactors::Tuple{Int, Int, Int}`: The factor you want to replicate the box by # Returns -- `box::Box`: Fully formed Box object + + - `box::Box`: Fully formed Box object """ function replicate(box::Box, repfactors::Tuple{Int, Int, Int}) - return Box(box.a * repfactors[1], box.b * repfactors[2], box.c * repfactors[3], - box.α, box.β, box.γ) + return Box( + box.a * repfactors[1], + box.b * repfactors[2], + box.c * repfactors[3], + box.α, + box.β, + box.γ + ) end """ @@ -161,21 +174,29 @@ same as the crystal structure filename but with a .vtk extension. Appends ".vtk" extension to `filename` automatically if not passed. # Arguments -- `box::Box`: a Bravais lattice -- `filename::AbstractString`: filename of the .vtk file output (absolute path) -- `crystal::Crystal`: a crystal structure object -- `center_at_origin::Bool`: center box at origin if true. if false, the origin is the corner of the box. -- `verbose::Bool`: if true, print the name of the written file to the console. + + - `box::Box`: a Bravais lattice + - `filename::AbstractString`: filename of the .vtk file output (absolute path) + - `crystal::Crystal`: a crystal structure object + - `center_at_origin::Bool`: center box at origin if true. if false, the origin is the corner of the box. + - `verbose::Bool`: if true, print the name of the written file to the console. """ -function write_vtk(box::Box, filename::AbstractString; verbose::Bool=false, - center_at_origin::Bool=false) - if ! occursin(".vtk", filename) +function write_vtk( + box::Box, + filename::AbstractString; + verbose::Bool=false, + center_at_origin::Bool=false +) + if !occursin(".vtk", filename) filename *= ".vtk" end open(filename, "w") do vtk_file - @printf(vtk_file, "# vtk DataFile Version 2.0\nunit cell boundary\n - ASCII\nDATASET POLYDATA\nPOINTS 8 double\n") + @printf( + vtk_file, + "# vtk DataFile Version 2.0\nunit cell boundary\n + ASCII\nDATASET POLYDATA\nPOINTS 8 double\n" + ) x_shift = zeros(3) if center_at_origin @@ -183,26 +204,33 @@ function write_vtk(box::Box, filename::AbstractString; verbose::Bool=false, end # write points on boundary of unit cell - for i = 0:1 - for j = 0:1 - for k = 0:1 + for i in 0:1 + for j in 0:1 + for k in 0:1 xf = [i, j, k] # fractional coordinates of corner cornerpoint = box.f_to_c * xf - x_shift - @printf(vtk_file, "%.3f %.3f %.3f\n", - cornerpoint[1], cornerpoint[2], cornerpoint[3]) + @printf( + vtk_file, + "%.3f %.3f %.3f\n", + cornerpoint[1], + cornerpoint[2], + cornerpoint[3] + ) end end end # define connections - @printf(vtk_file, "LINES 12 36\n2 0 1\n2 0 2\n2 1 3\n2 2 3\n2 4 5\n2 4 6\n2 5 7\n2 6 7\n2 0 4\n2 1 5\n2 2 6\n2 3 7\n") + @printf( + vtk_file, + "LINES 12 36\n2 0 1\n2 0 2\n2 1 3\n2 2 3\n2 4 5\n2 4 6\n2 5 7\n2 6 7\n2 0 4\n2 1 5\n2 2 6\n2 3 7\n" + ) end if verbose println("See ", filename) end end - """ inside_box = inside(frac_coords) # true or false inside_box = inside(cart_coords, box) @@ -215,10 +243,11 @@ returns true if coords are all inside a box and false otherwise. if a molecule or crystal is passed, both atoms and charges must be inside the box. # Arguments -* `coords::Coords` the coordinates (works for `Cart` and `Frac`) -* `molecule::Molecule{T}`: a molecule in either `T::Frac` or `T::Cart` coords -* `crystal::Crystal`: a crystal -* `box::Box` the box (only needed if Cartesian) + + - `coords::Coords` the coordinates (works for `Cart` and `Frac`) + - `molecule::Molecule{T}`: a molecule in either `T::Frac` or `T::Cart` coords + - `crystal::Crystal`: a crystal + - `box::Box` the box (only needed if Cartesian) """ inside(coords::Frac) = all(coords.xf .<= 1.0) && all(coords.xf .>= 0.0) inside(coords::Cart, box::Box) = inside(Frac(coords, box)) @@ -229,22 +258,39 @@ translate_by!(coords::Cart, dx::Frac, box::Box) = translate_by!(coords, Cart(dx, function Base.show(io::IO, box::Box) println(io, "Bravais unit cell of a crystal.") - @printf(io, "\tUnit cell angles α = %f deg. β = %f deg. γ = %f deg.\n", - box.α * 180.0 / π, box.β * 180.0 / π, box.γ * 180.0 / π) - @printf(io, "\tUnit cell dimensions a = %f Å. b = %f Å, c = %f Å\n", - box.a, box.b, box.c) + @printf( + io, + "\tUnit cell angles α = %f deg. β = %f deg. γ = %f deg.\n", + box.α * 180.0 / π, + box.β * 180.0 / π, + box.γ * 180.0 / π + ) + @printf( + io, + "\tUnit cell dimensions a = %f Å. b = %f Å, c = %f Å\n", + box.a, + box.b, + box.c + ) @printf(io, "\tVolume of unit cell: %f ų\n", box.Ω) end -function Base.isapprox(box1::Box, box2::Box; atol::Real=0.0, rtol::Real=atol > 0 ? 0.0 : sqrt(eps())) - return (isapprox(box1.a, box2.a, atol=atol, rtol=rtol) && - isapprox(box1.b, box2.b, atol=atol, rtol=rtol) && - isapprox(box1.c, box2.c, atol=atol, rtol=rtol) && - isapprox(box1.α, box2.α, atol=atol, rtol=rtol) && - isapprox(box1.β, box2.β, atol=atol, rtol=rtol) && - isapprox(box1.γ, box2.γ, atol=atol, rtol=rtol) && - isapprox(box1.Ω, box2.Ω, atol=atol, rtol=rtol) && - isapprox(box1.f_to_c, box2.f_to_c, atol=atol, rtol=rtol) && - isapprox(box1.c_to_f, box2.c_to_f, atol=atol, rtol=rtol) && - isapprox(box1.reciprocal_lattice, box2.reciprocal_lattice, atol=atol, rtol=rtol)) +function Base.isapprox( + box1::Box, + box2::Box; + atol::Real=0.0, + rtol::Real=atol > 0 ? 0.0 : sqrt(eps()) +) + return ( + isapprox(box1.a, box2.a; atol=atol, rtol=rtol) && + isapprox(box1.b, box2.b; atol=atol, rtol=rtol) && + isapprox(box1.c, box2.c; atol=atol, rtol=rtol) && + isapprox(box1.α, box2.α; atol=atol, rtol=rtol) && + isapprox(box1.β, box2.β; atol=atol, rtol=rtol) && + isapprox(box1.γ, box2.γ; atol=atol, rtol=rtol) && + isapprox(box1.Ω, box2.Ω; atol=atol, rtol=rtol) && + isapprox(box1.f_to_c, box2.f_to_c; atol=atol, rtol=rtol) && + isapprox(box1.c_to_f, box2.c_to_f; atol=atol, rtol=rtol) && + isapprox(box1.reciprocal_lattice, box2.reciprocal_lattice; atol=atol, rtol=rtol) + ) end diff --git a/src/crystal.jl b/src/crystal.jl index d1f5118..f504f78 100644 --- a/src/crystal.jl +++ b/src/crystal.jl @@ -2,13 +2,14 @@ SymmetryInfo(symmetry, space_group, is_p1) # Attributes -- `operations::Array{Function, 2}`: 2D array of anonymous functions that represent + + - `operations::Array{Function, 2}`: 2D array of anonymous functions that represent the symmetry operations. If the structure is in P1 there will be one symmetry operation. -- `space_group::AbstractString`: The name of the space group. This is stored + - `space_group::AbstractString`: The name of the space group. This is stored so that it can be written out again in the write_cif function. The space group is not used to verify the symmetry rules. -- `is_p1::Bool`: Stores whether the crystal is currently in P1 symmetry. + - `is_p1::Bool`: Stores whether the crystal is currently in P1 symmetry. """ struct SymmetryInfo operations::Array{String, 2} @@ -27,8 +28,9 @@ struct Crystal <: AbstractSystem{3} end # default constructor without bond info or symmetry info -Crystal(name::String, box::Box, atoms::Atoms{Frac}, charges::Charges{Frac}) = Crystal( - name, box, atoms, charges, MetaGraph(atoms.n), SymmetryInfo()) +function Crystal(name::String, box::Box, atoms::Atoms{Frac}, charges::Charges{Frac}) + return Crystal(name, box, atoms, charges, MetaGraph(atoms.n), SymmetryInfo()) +end const NET_CHARGE_TOL = 1e-4 # net charge tolerance @@ -48,53 +50,58 @@ Read a crystal structure file (.cif or .cssr) and populate a `Crystal` data stru or construct a `Crystal` data structure directly. # Arguments -- `filename::String`: the name of the crystal structure file (include ".cif" or ".cssr") read from `PATH_TO_CRYSTALS`. -- `check_neutrality::Bool`: check for charge neutrality -- `net_charge_tol::Float64`: when checking for charge neutrality, throw an error if the absolute value of the net charge is larger than this value. -- `check_overlap::Bool`: throw an error if overlapping atoms are detected. -- `convert_to_p1::Bool`: If the structure is not in P1 it will be converted to + + - `filename::String`: the name of the crystal structure file (include ".cif" or ".cssr") read from `PATH_TO_CRYSTALS`. + - `check_neutrality::Bool`: check for charge neutrality + - `net_charge_tol::Float64`: when checking for charge neutrality, throw an error if the absolute value of the net charge is larger than this value. + - `check_overlap::Bool`: throw an error if overlapping atoms are detected. + - `convert_to_p1::Bool`: If the structure is not in P1 it will be converted to P1 symmetry using the symmetry rules from the `_symmetry_equiv_pos_as_xyz` list in the .cif file. (We do not use the space groups name to look up symmetry rules). -- `read_bonds_from_file::Bool`: Whether or not to read bonding information from + - `read_bonds_from_file::Bool`: Whether or not to read bonding information from cif file. If false, the bonds can be inferred later. note that, if the crystal is not in P1 symmetry, we cannot *both* read bonds and convert to P1 symmetry. -- `wrap_coords::Bool`: if `true`, enforce that fractional coords of atoms and charges are in [0,1]³ by mod(x, 1) -- `include_zero_charges::Bool`: if `false`, do not include in `crystal.charges` atoms which have zero charges, in order to speed up the electrostatic calculations. + - `wrap_coords::Bool`: if `true`, enforce that fractional coords of atoms and charges are in [0,1]³ by mod(x, 1) + - `include_zero_charges::Bool`: if `false`, do not include in `crystal.charges` atoms which have zero charges, in order to speed up the electrostatic calculations. If `true,` include the atoms in `crystal.charges` that have zero charge, ensuring that the number of atoms is equal to the number of charges and that `crystal.charges.coords.xf` and `crystal.atoms.coords.xf` are the same. -- `remove_duplicates::Bool`: remove duplicate atoms and charges. an atom is duplicate only if it is the same species. -- `species_col::Array{String}`: which column to use for species identification for `crystal.atoms.species`. we use a priority list: + - `remove_duplicates::Bool`: remove duplicate atoms and charges. an atom is duplicate only if it is the same species. + - `species_col::Array{String}`: which column to use for species identification for `crystal.atoms.species`. we use a priority list: we check for the first entry of `species_col` in the .cif file; if not present, we then use the second entry, and so on. -- `infer_bonds::Bool`: if true, bonds are inferred automatically. If set, must specify `periodic_boundaries`. By default, bonds are not inferred. -- `periodic_boundaries::Union{Bool, Missing}`: use with `infer_bonds` to specify treatment of the unit cell boundary. Set `true` to treat the unit cell edge as a periodic boundary (allow bonds across it); set `false` to restrict bonding to within the local unit cell. By default, no bonds are inferred. -- `absolute_path::Bool`: set `true` to load from an absolute path instead of the default path stored in `rc[:paths][:crystals]`. + - `infer_bonds::Bool`: if true, bonds are inferred automatically. If set, must specify `periodic_boundaries`. By default, bonds are not inferred. + - `periodic_boundaries::Union{Bool, Missing}`: use with `infer_bonds` to specify treatment of the unit cell boundary. Set `true` to treat the unit cell edge as a periodic boundary (allow bonds across it); set `false` to restrict bonding to within the local unit cell. By default, no bonds are inferred. + - `absolute_path::Bool`: set `true` to load from an absolute path instead of the default path stored in `rc[:paths][:crystals]`. # Returns -- `crystal::Crystal`: A crystal containing the crystal structure information + + - `crystal::Crystal`: A crystal containing the crystal structure information # Attributes -- `name::AbstractString`: name of crystal structure -- `box::Box`: unit cell (Bravais Lattice) -- `atoms::Atoms`: list of Atoms in crystal unit cell -- `charges::Charges`: list of point charges in crystal unit cell -- `bonds::MetaGraph`: Unweighted, undirected graph showing all of the atoms + + - `name::AbstractString`: name of crystal structure + - `box::Box`: unit cell (Bravais Lattice) + - `atoms::Atoms`: list of Atoms in crystal unit cell + - `charges::Charges`: list of point charges in crystal unit cell + - `bonds::MetaGraph`: Unweighted, undirected graph showing all of the atoms that are bonded within the crystal -- `symmetry::SymmetryInfo`: symmetry inforomation + - `symmetry::SymmetryInfo`: symmetry inforomation """ -function Crystal(filename::String; - check_neutrality::Bool=true, - net_charge_tol::Float64=NET_CHARGE_TOL, - check_overlap::Bool=true, - convert_to_p1::Bool=true, - read_bonds_from_file::Bool=false, - wrap_coords::Bool=true, - include_zero_charges::Bool=true, - remove_duplicates::Bool=false, - species_col::Array{String, 1}=["_atom_site_type_symbol", "_atom_site_label"], - infer_bonds::Bool=false, - periodic_boundaries::Union{Bool, Missing}=missing, - absolute_path::Bool=false) +function Crystal( + filename::String; + check_neutrality::Bool=true, + net_charge_tol::Float64=NET_CHARGE_TOL, + check_overlap::Bool=true, + convert_to_p1::Bool=true, + read_bonds_from_file::Bool=false, + wrap_coords::Bool=true, + include_zero_charges::Bool=true, + remove_duplicates::Bool=false, + species_col::Array{String, 1}=["_atom_site_type_symbol", "_atom_site_label"], + infer_bonds::Bool=false, + periodic_boundaries::Union{Bool, Missing}=missing, + absolute_path::Bool=false +) # Read file extension. Ensure we can read the file type extension = split(filename, ".")[end] - if ! (extension in ["cif", "cssr"]) + if !(extension in ["cif", "cssr"]) error("I can only read .cif or .cssr crystal structure files.") end @@ -123,7 +130,6 @@ function Crystal(filename::String; is_p1 = false space_group = "" - # Start of .cif reader ################################### # CIF READER @@ -150,13 +156,13 @@ function Crystal(filename::String; end # Make sure the space group is P1 - if line[1] == "_symmetry_space_group_name_H-M" || line[1] == "_space_group_name_Hall" + if line[1] == "_symmetry_space_group_name_H-M" || + line[1] == "_space_group_name_Hall" # use anonymous function to combine all terms past the first # to extract space group name space_group = reduce((x, y) -> x * " " * y, line[2:end]) - space_group = split(space_group, [''', '"'], keepempty=false)[1] - if space_group == "P1" || space_group == "P 1" || - space_group == "-P1" + space_group = split(space_group, [''', '"']; keepempty=false)[1] + if space_group == "P1" || space_group == "P 1" || space_group == "-P1" # simplify by only having one P1 space_group name space_group = "P1" is_p1 = true @@ -179,19 +185,25 @@ function Crystal(filename::String; i += 1 end - fractional = fractional || haskey(name_to_column, "_atom_site_fract_x") && - haskey(name_to_column, "_atom_site_fract_y") && - haskey(name_to_column, "_atom_site_fract_z") + fractional = + fractional || + haskey(name_to_column, "_atom_site_fract_x") && + haskey(name_to_column, "_atom_site_fract_y") && + haskey(name_to_column, "_atom_site_fract_z") # if the file provides cartesian coordinates - cartesian = cartesian || ! fractional && haskey(name_to_column, "_atom_site_Cartn_x") && - haskey(name_to_column, "_atom_site_Cartn_y") && - haskey(name_to_column, "_atom_site_Cartn_z") - # if both are provided, will default - # to using fractional, so keep cartesian - # false + cartesian = + cartesian || + !fractional && + haskey(name_to_column, "_atom_site_Cartn_x") && + haskey(name_to_column, "_atom_site_Cartn_y") && + haskey(name_to_column, "_atom_site_Cartn_z") + # if both are provided, will default + # to using fractional, so keep cartesian + # false # Assign species_column by matching to priority list - if haskey(name_to_column, "_atom_site_Cartn_x") || haskey(name_to_column, "_atom_site_fract_x") # to have entered _atom_site loop + if haskey(name_to_column, "_atom_site_Cartn_x") || + haskey(name_to_column, "_atom_site_fract_x") # to have entered _atom_site loop found_species_col = false for col in species_col if col ∈ keys(name_to_column) @@ -200,7 +212,7 @@ function Crystal(filename::String; break end end - if ! found_species_col + if !found_species_col @error "could not find species_col=$(species_col) columns in the .cif file for atomic species labels for $(filename)." end end @@ -217,9 +229,12 @@ function Crystal(filename::String; # not really one column) the length(name_to_column) + 2 # should catch this hopefully there aren't other weird # ways of writing cifs... - while i <= length(lines) && length(lines[i]) > 0 && lines[i][1] != '_' && !occursin("loop_", lines[i]) + while i <= length(lines) && + length(lines[i]) > 0 && + lines[i][1] != '_' && + !occursin("loop_", lines[i]) line = lines[i] - sym_funcs = split(line, [' ', ',', ''', '"', '\t'], keepempty=false) + sym_funcs = split(line, [' ', ',', ''', '"', '\t']; keepempty=false) if length(sym_funcs) == 0 break end @@ -229,7 +244,7 @@ function Crystal(filename::String; new_sym_rule = Array{AbstractString, 1}(undef, 3) sym_start = name_to_column["_symmetry_equiv_pos_as_xyz"] - 1 - for j = 1:3 + for j in 1:3 new_sym_rule[j] = sym_funcs[j + sym_start] end @@ -243,28 +258,53 @@ function Crystal(filename::String; # finish reading in symmetry information, skip to next # iteration of outer while-loop continue - # ===================== - # FRACTIONAL READER - # ===================== - elseif fractional && ! atom_info + # ===================== + # FRACTIONAL READER + # ===================== + elseif fractional && !atom_info atom_info = true - while i <= length(lines) && length(split(lines[i])) == length(name_to_column) + while i <= length(lines) && + length(split(lines[i])) == length(name_to_column) line = split(lines[i]) push!(species, Symbol(line[name_to_column[species_column]])) - coords = [coords [mod(parse(Float64, split(line[name_to_column["_atom_site_fract_x"]], '(')[1]), 1.0), - mod(parse(Float64, split(line[name_to_column["_atom_site_fract_y"]], '(')[1]), 1.0), - mod(parse(Float64, split(line[name_to_column["_atom_site_fract_z"]], '(')[1]), 1.0)]] + coords = [coords [ + mod( + parse( + Float64, + split(line[name_to_column["_atom_site_fract_x"]], '(')[1] + ), + 1.0 + ), + mod( + parse( + Float64, + split(line[name_to_column["_atom_site_fract_y"]], '(')[1] + ), + 1.0 + ), + mod( + parse( + Float64, + split(line[name_to_column["_atom_site_fract_z"]], '(')[1] + ), + 1.0 + ) + ]] # If charges present, import them if haskey(name_to_column, "_atom_site_charge") - push!(charge_values, parse(Float64, line[name_to_column["_atom_site_charge"]])) + push!( + charge_values, + parse(Float64, line[name_to_column["_atom_site_charge"]]) + ) else push!(charge_values, NaN) end # add to label_num_to_idx so that bonds can be converted later if read_bonds_from_file - label_num_to_idx[line[name_to_column[species_column]]] = length(species) + label_num_to_idx[line[name_to_column[species_column]]] = + length(species) end # iterate to next line in file i += 1 @@ -275,28 +315,44 @@ function Crystal(filename::String; # iteration of outer while-loop # prevents skipping a line after finishing reading atoms continue - # ===================== - # CARTESIAN READER - # ===================== - elseif cartesian && ! atom_info + # ===================== + # CARTESIAN READER + # ===================== + elseif cartesian && !atom_info atom_info = true - while i <= length(lines) && length(split(lines[i])) == length(name_to_column) + while i <= length(lines) && + length(split(lines[i])) == length(name_to_column) line = split(lines[i]) push!(species, Symbol(line[name_to_column[species_column]])) - coords = [coords [parse(Float64, split(line[name_to_column["_atom_site_Cartn_x"]], '(')[1]), - parse(Float64, split(line[name_to_column["_atom_site_Cartn_y"]], '(')[1]), - parse(Float64, split(line[name_to_column["_atom_site_Cartn_z"]], '(')[1])]] + coords = [coords [ + parse( + Float64, + split(line[name_to_column["_atom_site_Cartn_x"]], '(')[1] + ), + parse( + Float64, + split(line[name_to_column["_atom_site_Cartn_y"]], '(')[1] + ), + parse( + Float64, + split(line[name_to_column["_atom_site_Cartn_z"]], '(')[1] + ) + ]] # If charges present, import them if haskey(name_to_column, "_atom_site_charge") - push!(charge_values, parse(Float64, line[name_to_column["_atom_site_charge"]])) + push!( + charge_values, + parse(Float64, line[name_to_column["_atom_site_charge"]]) + ) else push!(charge_values, NaN) end # add to label_num_to_idx so that bonds can be converted later if read_bonds_from_file - label_num_to_idx[line[name_to_column[species_column]]] = length(species) + label_num_to_idx[line[name_to_column[species_column]]] = + length(species) end # iterate to next line in file i += 1 @@ -307,18 +363,26 @@ function Crystal(filename::String; # iteration of outer while-loop # prevents skipping a line after finishing reading atoms continue - # ===================== - # BOND READER - # ===================== + # ===================== + # BOND READER + # ===================== elseif read_bonds_from_file && haskey(name_to_column, "_geom_bond_atom_site_label_1") && haskey(name_to_column, "_geom_bond_atom_site_label_2") - while i <= length(lines) && length(split(lines[i])) == length(name_to_column) + while i <= length(lines) && + length(split(lines[i])) == length(name_to_column) line = split(lines[i]) - atom_one_idx = label_num_to_idx[line[name_to_column["_geom_bond_atom_site_label_1"]]] - atom_two_idx = label_num_to_idx[line[name_to_column["_geom_bond_atom_site_label_2"]]] - add_edge!(bonds, atom_one_idx, atom_two_idx, :distance, - distance(coords, box, atom_one_idx, atom_two_idx, true)) + atom_one_idx = + label_num_to_idx[line[name_to_column["_geom_bond_atom_site_label_1"]]] + atom_two_idx = + label_num_to_idx[line[name_to_column["_geom_bond_atom_site_label_2"]]] + add_edge!( + bonds, + atom_one_idx, + atom_two_idx, + :distance, + distance(coords, box, atom_one_idx, atom_two_idx, true) + ) # iterate to next line in file i += 1 end @@ -331,14 +395,14 @@ function Crystal(filename::String; # pick up unit cell lengths for axis in ["a", "b", "c"] if line[1] == @sprintf("_cell_length_%s", axis) - data[axis] = parse(Float64, split(line[2],'(')[1]) + data[axis] = parse(Float64, split(line[2], '(')[1]) end end # pick up unit cell angles for angle in ["alpha", "beta", "gamma"] if line[1] == @sprintf("_cell_angle_%s", angle) - data[angle] = parse(Float64, split(line[2],'(')[1]) * pi / 180.0 + data[angle] = parse(Float64, split(line[2], '(')[1]) * pi / 180.0 end end @@ -350,9 +414,14 @@ function Crystal(filename::String; end # Structure must either be in P1 symmetry or have replication information - if ! is_p1 && !symmetry_info - error(@sprintf("%s is not in P1 symmetry and the .cif does not have a _symmetry_equiv_pos_as_xyz column - for us to apply symmetry operations to convert into P1 symmetry.", filename)) + if !is_p1 && !symmetry_info + error( + @sprintf( + "%s is not in P1 symmetry and the .cif does not have a _symmetry_equiv_pos_as_xyz column + for us to apply symmetry operations to convert into P1 symmetry.", + filename + ) + ) end a = data["a"] @@ -363,14 +432,14 @@ function Crystal(filename::String; γ = data["gamma"] # redo coordinates if they were read in cartesian - if cartesian && ! fractional + if cartesian && !fractional coords = Box(a, b, c, α, β, γ).c_to_f * coords end - # Start of cssr reader #TODO make sure this works for different .cssr files! - ################################### - # CSSR READER - ################################### + # Start of cssr reader #TODO make sure this works for different .cssr files! + ################################### + # CSSR READER + ################################### elseif extension == "cssr" # First line contains unit cell lenghts line = split(lines[1]) @@ -388,7 +457,7 @@ function Crystal(filename::String; bonds = MetaGraph(n_atoms) # Read in atoms and fractional coordinates - for i = 1:n_atoms + for i in 1:n_atoms line = split(lines[4 + i]) push!(species, Symbol(line[2])) @@ -399,8 +468,8 @@ function Crystal(filename::String; push!(charge_values, parse(Float64, line[14])) end - for i = 1:n_atoms - coords = [ coords [xf[i], yf[i], zf[i]] ] + for i in 1:n_atoms + coords = [coords [xf[i], yf[i], zf[i]]] end # add P1 symmetry rules for consistency @@ -414,10 +483,10 @@ function Crystal(filename::String; # construct atoms attribute of crystal atoms = Atoms(species, Frac(coords)) # construct charges attribute of crystal - if all(isnan.(charge_values)) + if all(isnan.(charge_values)) # if no charges column, charge_values will be all NaN charges = Charges{Frac}(0) - elseif ! include_zero_charges + elseif !include_zero_charges # include only nonzero charges idx_nz = charge_values .!= 0.0 charges = Charges(charge_values[idx_nz], Frac(coords[:, idx_nz])) @@ -429,10 +498,13 @@ function Crystal(filename::String; symmetry = SymmetryInfo(operations, space_group, is_p1) crystal = Crystal(filename, box, atoms, charges, bonds, symmetry) - if convert_to_p1 && ! is_p1 && ! read_bonds_from_file - @info @sprintf("Crystal %s has %s space group. I am converting it to P1 symmetry. - To prevent this, pass `convert_to_p1=false` to the `Crystal` constructor.\n", - crystal.name, crystal.symmetry.space_group) + if convert_to_p1 && !is_p1 && !read_bonds_from_file + @info @sprintf( + "Crystal %s has %s space group. I am converting it to P1 symmetry. +To prevent this, pass `convert_to_p1=false` to the `Crystal` constructor.\n", + crystal.name, + crystal.symmetry.space_group + ) crystal = apply_symmetry_operations(crystal) end @@ -445,10 +517,16 @@ function Crystal(filename::String; end if check_neutrality - if ! neutral(crystal, net_charge_tol) - error(@sprintf("Crystal %s is not charge neutral; net charge is %f e. Ignore - this error message by passing check_neutrality=false or increasing the - net charge tolerance `net_charge_tol`\n", crystal.name, net_charge(crystal))) + if !neutral(crystal, net_charge_tol) + error( + @sprintf( + "Crystal %s is not charge neutral; net charge is %f e. Ignore + this error message by passing check_neutrality=false or increasing the + net charge tolerance `net_charge_tol`\n", + crystal.name, + net_charge(crystal) + ) + ) end end @@ -475,7 +553,6 @@ function Crystal(filename::String; return crystal end - """ renamed_xtal = rename_xtal(xtal, "new name") @@ -486,18 +563,17 @@ function rename_xtal(xtal::Crystal, name::String)::Crystal nothing_xtal = Crystal( "nothing", xtal.box, - Atoms(0, Symbol[], Frac(zeros(3,0))), + Atoms(0, Symbol[], Frac(zeros(3, 0))), Charges{Frac}(0) ) # add the "nothing" crystal to the input crystal and change the name - return +(xtal, nothing_xtal, name=name) + return +(xtal, nothing_xtal; name=name) end - # documented in matter.jl function wrap!(crystal::Crystal) wrap!(crystal.atoms.coords) - wrap!(crystal.charges.coords) + return wrap!(crystal.charges.coords) end # documented in matter.jl @@ -511,11 +587,13 @@ replicate the atoms and charges in a `Crystal` in positive directions to constru coordinates will be rescaled to be in [0, 1]. # arguments -- `crystal::Crystal`: The crystal to replicate -- `repfactors::Tuple{Int, Int, Int}`: The factors by which to replicate the crystal structure in each crystallographic direction (a, b, c). + + - `crystal::Crystal`: The crystal to replicate + - `repfactors::Tuple{Int, Int, Int}`: The factors by which to replicate the crystal structure in each crystallographic direction (a, b, c). # returns -- `replicated_frame::Crystal`: replicated crystal + + - `replicated_frame::Crystal`: replicated crystal """ function replicate(crystal::Crystal, repfactors::Tuple{Int, Int, Int}) if ne(crystal.bonds) != 0 @@ -535,11 +613,14 @@ function replicate(crystal::Crystal, repfactors::Tuple{Int, Int, Int}) atom_counter = 0 charge_counter = 0 - for ra = 0:(repfactors[1] - 1), rb = 0:(repfactors[2] - 1), rc = 0:(repfactors[3] - 1) + for ra in 0:(repfactors[1] - 1), + rb in 0:(repfactors[2] - 1), + rc in 0:(repfactors[3] - 1) + xf_shift = 1.0 * [ra, rb, rc] # replicate atoms - for i = 1:crystal.atoms.n + for i in 1:(crystal.atoms.n) atom_counter += 1 atoms.species[atom_counter] = crystal.atoms.species[i] @@ -549,7 +630,7 @@ function replicate(crystal::Crystal, repfactors::Tuple{Int, Int, Int}) end # replicate charges - for i = 1:crystal.charges.n + for i in 1:(crystal.charges.n) charge_counter += 1 charges.q[charge_counter] = crystal.charges.q[i] @@ -563,27 +644,32 @@ function replicate(crystal::Crystal, repfactors::Tuple{Int, Int, Int}) end # doc string in Misc.jl -xyz_filename(crystal::Crystal) = replace(replace(crystal.name, ".cif" => ""), ".cssr" => "") * ".xyz" +function xyz_filename(crystal::Crystal) + return replace(replace(crystal.name, ".cif" => ""), ".cssr" => "") * ".xyz" +end -function write_xyz(crystal::Crystal, filename::AbstractString=xyz_filename(crystal); center_at_origin::Bool=false, kwargs...) - atoms = Atoms(crystal.atoms.species, - Cart(crystal.atoms.coords, crystal.box) - ) # put in Cartesian +function write_xyz( + crystal::Crystal, + filename::AbstractString=xyz_filename(crystal); + center_at_origin::Bool=false, + kwargs... +) + atoms = Atoms(crystal.atoms.species, Cart(crystal.atoms.coords, crystal.box)) # put in Cartesian if center_at_origin x_c = crystal.box.f_to_c * [0.5, 0.5, 0.5] atoms.coords.x .-= x_c - end - write_xyz(atoms, filename; kwargs...) + return write_xyz(atoms, filename; kwargs...) end -write_vtk(crystal::Crystal, filename::AbstractString; kwargs...) = write_vtk(crystal.box, filename; kwargs...) +function write_vtk(crystal::Crystal, filename::AbstractString; kwargs...) + return write_vtk(crystal.box, filename; kwargs...) +end write_vtk(crystal::Crystal; kwargs...) = write_vtk(crystal.box, crystal.name; kwargs...) # docstring in matter.jl neutral(crystal::Crystal, tol::Float64=1e-5) = neutral(crystal.charges, tol) - """ charged = has_charges(crystal) # true or false charged = has_charges(molecule) # true or false @@ -592,7 +678,6 @@ neutral(crystal::Crystal, tol::Float64=1e-5) = neutral(crystal.charges, tol) """ has_charges(crystal::Crystal) = crystal.charges.n > 0 - # e.g. `:Ca23` -> `:Ca` function strip_number_from_label(atom_label::Symbol) atom_label_string = String(atom_label) @@ -602,11 +687,10 @@ function strip_number_from_label(atom_label::Symbol) # nothing to strip return atom_label else - return Symbol(atom_label_string[1:findfirst(.! isletter_vector) - 1]) + return Symbol(atom_label_string[1:(findfirst(.!isletter_vector) - 1)]) end end - """ strip_numbers_from_atom_labels!(crystal) @@ -614,18 +698,21 @@ Strip numbers from labels for `crystal.atoms`. Precisely, for `atom` in `crystal.atoms`, find the first number that appears in `atom`. Remove this number and all following characters from `atom`. e.g. C12 --> C - Ba12A_3 --> Ba +Ba12A_3 --> Ba # Arguments -- `crystal::Crystal`: The crystal containing the crystal structure information + + - `crystal::Crystal`: The crystal containing the crystal structure information """ function strip_numbers_from_atom_labels!(crystal::Crystal) - for i = 1:crystal.atoms.n + for i in 1:(crystal.atoms.n) crystal.atoms.species[i] = strip_number_from_label(crystal.atoms.species[i]) end end -vtk_filename(crystal::Crystal) = replace(replace(crystal.name, ".cif" => ""), ".cssr" => "") * ".vtk" +function vtk_filename(crystal::Crystal) + return replace(replace(crystal.name, ".cif" => ""), ".cssr" => "") * ".vtk" +end """ formula = empirical_formula(crystal, verbose=false) @@ -633,17 +720,19 @@ vtk_filename(crystal::Crystal) = replace(replace(crystal.name, ".cif" => ""), ". Find the irreducible chemical formula of a crystal structure. # Arguments -- `crystal::Crystal`: The crystal containing the crystal structure information -- `verbose::Bool`: If `true`, will print the chemical formula as well + + - `crystal::Crystal`: The crystal containing the crystal structure information + - `verbose::Bool`: If `true`, will print the chemical formula as well # Returns -- `formula::Dict{Symbol, Int}`: A dictionary with the irreducible chemical formula of a crystal structure + + - `formula::Dict{Symbol, Int}`: A dictionary with the irreducible chemical formula of a crystal structure """ function empirical_formula(crystal::Crystal; verbose::Bool=false) unique_atoms = unique(crystal.atoms.species) # use dictionary to count atom types atom_counts = Dict{Symbol, Int}([a => 0 for a in unique_atoms]) - for i = 1:crystal.atoms.n + for i in 1:(crystal.atoms.n) atom_counts[crystal.atoms.species[i]] += 1 end @@ -668,22 +757,23 @@ function empirical_formula(crystal::Crystal; verbose::Bool=false) end """ - mass_of_crystal = molecular_weight(crystal) Calculates the molecular weight of a unit cell of the crystal in amu using information stored in `data/atomicmasses.csv`. # Arguments -- `crystal::Crystal`: The crystal containing the crystal structure information + + - `crystal::Crystal`: The crystal containing the crystal structure information # Returns -- `mass_of_crystal::Float64`: The molecular weight of a unit cell of the crystal in amu + + - `mass_of_crystal::Float64`: The molecular weight of a unit cell of the crystal in amu """ function molecular_weight(crystal::Crystal) atomic_masses = rc[:atomic_masses] mass = 0.0 - for i = 1:crystal.atoms.n + for i in 1:(crystal.atoms.n) mass += atomic_masses[crystal.atoms.species[i]] end @@ -696,10 +786,12 @@ end Compute the crystal density of a crystal. Pulls atomic masses from [`read_atomic_masses`](@ref). # Arguments -- `crystal::Crystal`: The crystal containing the crystal structure information + + - `crystal::Crystal`: The crystal containing the crystal structure information # Returns -- `ρ::Float64`: The crystal density of a crystal in kg/m³ + + - `ρ::Float64`: The crystal density of a crystal in kg/m³ # --> kg/m3 """ crystal_density(crystal::Crystal) = molecular_weight(crystal) / crystal.box.Ω * 1660.53892 # --> kg/m3 @@ -709,10 +801,12 @@ crystal_density(crystal::Crystal) = molecular_weight(crystal) / crystal.box.Ω * Convert a crystal to P1 symmetry based on internal symmetry rules. This will return a new crystal. # Arguments -- `non_p1_crystal::Crystal`: The crystal to be converted to P1 symmetry + + - `non_p1_crystal::Crystal`: The crystal to be converted to P1 symmetry # Returns -- `P1_crystal::Crystal`: The crystal after it has been converted to P1 + + - `P1_crystal::Crystal`: The crystal after it has been converted to P1 symmetry. The new symmetry rules will be the P1 symmetry rules """ function apply_symmetry_operations(crystal::Crystal) @@ -732,42 +826,44 @@ function apply_symmetry_operations(crystal::Crystal) for sr in 1:nb_symmetry_ops # loop over all atoms in lower level symmetry sym_rule = eval.(Meta.parse.("(x, y, z) -> " .* crystal.symmetry.operations[:, sr])) - for a in 1:crystal.atoms.n + for a in 1:(crystal.atoms.n) # apply current symmetry rule to current atom for x, y, and z coordinates atom_id = (sr - 1) * crystal.atoms.n + a atoms.species[atom_id] = crystal.atoms.species[a] - atoms.coords.xf[:, atom_id] .= [Base.invokelatest.( - sym_rule[k], crystal.atoms.coords.xf[:, a]...) for k in 1:3] + atoms.coords.xf[:, atom_id] .= [ + Base.invokelatest.(sym_rule[k], crystal.atoms.coords.xf[:, a]...) for + k in 1:3 + ] end # loop over all charges in lower level symmetry - for c in 1:crystal.charges.n + for c in 1:(crystal.charges.n) # apply current symmetry rule to current atom for x, y, and z coordinates charge_id = (sr - 1) * crystal.charges.n + c charges.q[charge_id] = crystal.charges.q[c] - charges.coords.xf[:, charge_id] .= [Base.invokelatest.( - sym_rule[k], crystal.charges.coords.xf[:, c]...) for k in 1:3] + charges.coords.xf[:, charge_id] .= [ + Base.invokelatest.(sym_rule[k], crystal.charges.coords.xf[:, c]...) for + k in 1:3 + ] end end return Crystal(crystal.name, crystal.box, atoms, charges) end - """ assert_P1_symmetry(crystal::Crystal) Throw an error if and only if the crystal is not in P1 symmetry. """ function assert_P1_symmetry(crystal::Crystal) - if ! crystal.symmetry.is_p1 + if !crystal.symmetry.is_p1 error("the crystal " * crystal.name * " is not in P1 symmetry.\n To convert to P1 symmetry, try:\n \tcrystal_p1 = apply_symmetry_operations(crystal)") end end - """ write_cif(crystal, filename; fractional_coords=true, number_atoms=true) write_cif(crystal) # writes to file crystal.name[.cif] @@ -775,17 +871,23 @@ end Write a `crystal::Crystal` to a .cif file. # arguments -* `crystal::Crystal`: crystal to write to file -* `filename::String`: the filename of the `.cif` file. if ".cif" is not included as an extension, it will automatically be appended to the `filename` string. -* `fractional_coords::Bool=true`: write the coordinates of the atoms as fractional coords if `true`. if `false`, write Cartesian coords. -* `number_atoms::Bool=true`: write the atoms as "C1", "C2", "C3", ..., "N1", "N2", ... etc. to give each atom a unique identifier + + - `crystal::Crystal`: crystal to write to file + - `filename::String`: the filename of the `.cif` file. if ".cif" is not included as an extension, it will automatically be appended to the `filename` string. + - `fractional_coords::Bool=true`: write the coordinates of the atoms as fractional coords if `true`. if `false`, write Cartesian coords. + - `number_atoms::Bool=true`: write the atoms as "C1", "C2", "C3", ..., "N1", "N2", ... etc. to give each atom a unique identifier """ -function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=true, number_atoms::Bool=true) +function write_cif( + crystal::Crystal, + filename::String; + fractional_coords::Bool=true, + number_atoms::Bool=true +) if has_charges(crystal) if crystal.atoms.n != crystal.charges.n error("write_cif assumes equal numbers of Charges and Atoms (or zero charges)") end - if ! isapprox(crystal.charges.coords, crystal.atoms.coords) + if !isapprox(crystal.charges.coords, crystal.atoms.coords) error("write_cif needs coords of atoms and charges to correspond.") end end @@ -794,13 +896,13 @@ function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=t # create dictionary for tracking label numbers label_numbers = Dict{Symbol, Int}() for atom in crystal.atoms.species - if ! haskey(label_numbers, atom) + if !haskey(label_numbers, atom) label_numbers[atom] = 1 end end # append ".cif" to filename if it doesn't already have the extension - if ! occursin(".cif", filename) + if !occursin(".cif", filename) filename *= ".cif" end open(filename, "w") do cif_file @@ -812,7 +914,11 @@ function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=t @printf(cif_file, "data_%s_PM\n", split(crystal.name, ".")[1]) end - @printf(cif_file, "_symmetry_space_group_name_H-M '%s'\n", crystal.symmetry.space_group) + @printf( + cif_file, + "_symmetry_space_group_name_H-M '%s'\n", + crystal.symmetry.space_group + ) @printf(cif_file, "_cell_length_a %f\n", crystal.box.a) @printf(cif_file, "_cell_length_b %f\n", crystal.box.b) @@ -834,9 +940,15 @@ function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=t @printf(cif_file, "loop_\n_atom_site_label\n_atom_site_type_symbol\n") if fractional_coords - @printf(cif_file, "_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n") + @printf( + cif_file, + "_atom_site_fract_x\n_atom_site_fract_y\n_atom_site_fract_z\n" + ) else - @printf(cif_file, "_atom_site_Cartn_x\n_atom_site_Cartn_y\n_atom_site_Cartn_z\n") + @printf( + cif_file, + "_atom_site_Cartn_x\n_atom_site_Cartn_y\n_atom_site_Cartn_z\n" + ) end high_precision_charges = false # if, for neutrality, need high-precision charges if has_charges(crystal) @@ -850,20 +962,29 @@ function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=t end idx_to_label = Array{AbstractString, 1}(undef, crystal.atoms.n) - for i = 1:crystal.atoms.n + for i in 1:(crystal.atoms.n) # print label and type symbol - @printf(cif_file, "%s %s ", string(crystal.atoms.species[i]) * - (number_atoms ? string(label_numbers[crystal.atoms.species[i]]) : ""), - crystal.atoms.species[i]) + @printf( + cif_file, + "%s %s ", + string(crystal.atoms.species[i]) * + (number_atoms ? string(label_numbers[crystal.atoms.species[i]]) : ""), + crystal.atoms.species[i] + ) # store label for this atom idx - idx_to_label[i] = string(crystal.atoms.species[i]) * - string(label_numbers[crystal.atoms.species[i]]) + idx_to_label[i] = + string(crystal.atoms.species[i]) * + string(label_numbers[crystal.atoms.species[i]]) # increment label label_numbers[crystal.atoms.species[i]] += 1 if fractional_coords @printf(cif_file, "%f %f %f", crystal.atoms.coords.xf[:, i]...) else - @printf(cif_file, "%f %f %f", (crystal.box.f_to_c * crystal.atoms.coords.xf[:, i])...) + @printf( + cif_file, + "%f %f %f", + (crystal.box.f_to_c * crystal.atoms.coords.xf[:, i])... + ) end if has_charges(crystal) if high_precision_charges @@ -878,24 +999,35 @@ function write_cif(crystal::Crystal, filename::String; fractional_coords::Bool=t # only print bond information if it is in the crystal if ne(crystal.bonds) > 0 - if ! number_atoms + if !number_atoms error("must label atoms with numbers to write bond information.\n") end # print column names for bond information - @printf(cif_file, "\nloop_\n_geom_bond_atom_site_label_1\n_geom_bond_atom_site_label_2\n_geom_bond_distance\n") + @printf( + cif_file, + "\nloop_\n_geom_bond_atom_site_label_1\n_geom_bond_atom_site_label_2\n_geom_bond_distance\n" + ) for edge in collect(edges(crystal.bonds)) - dxf = crystal.atoms.coords.xf[:, edge.src] - crystal.atoms.coords.xf[:, edge.dst] + dxf = + crystal.atoms.coords.xf[:, edge.src] - + crystal.atoms.coords.xf[:, edge.dst] nearest_image!(dxf) - @printf(cif_file, "%s %s %0.5f\n", idx_to_label[edge.src], idx_to_label[edge.dst], - norm(dxf)) + @printf( + cif_file, + "%s %s %0.5f\n", + idx_to_label[edge.src], + idx_to_label[edge.dst], + norm(dxf) + ) end end end end -write_cif(crystal::Crystal; kwargs...) = write_cif(crystal, String(split(crystal.name, ".cif")[1]); kwargs...) - +function write_cif(crystal::Crystal; kwargs...) + return write_cif(crystal, String(split(crystal.name, ".cif")[1]); kwargs...) +end """ write_cssr(xtal, "myxtal.cssr") @@ -904,25 +1036,39 @@ write_cif(crystal::Crystal; kwargs...) = write_cif(crystal, String(split(crystal write crystal structure to `.cssr` format. # arguments -* `xtal::Crystal`: crystal to write to file -* `filename::String`: filename to write to. default is to write `.cssr` file to `pwd()`. will append `.cssr` if absent from `filename`. -* `quiet::Bool` : (optional) set `true` to suppress console output + + - `xtal::Crystal`: crystal to write to file + - `filename::String`: filename to write to. default is to write `.cssr` file to `pwd()`. will append `.cssr` if absent from `filename`. + - `quiet::Bool` : (optional) set `true` to suppress console output """ function write_cssr(xtal::Crystal, filename::String; quiet::Bool=false) @assert xtal.symmetry.is_p1 "crystal must be in P1 symmetry" - if ! occursin(".cssr", filename) + if !occursin(".cssr", filename) filename *= ".cssr" end open(filename, "w") do f @printf(f, "%f %f %f\n", xtal.box.a, xtal.box.b, xtal.box.c) - @printf(f, "%f %f %f SPGR = 1 P 1 OPT = 1\n", rad2deg(xtal.box.α), rad2deg(xtal.box.β), rad2deg(xtal.box.γ)) + @printf( + f, + "%f %f %f SPGR = 1 P 1 OPT = 1\n", + rad2deg(xtal.box.α), + rad2deg(xtal.box.β), + rad2deg(xtal.box.γ) + ) @printf(f, "%d 0\n0 xtal : xtal", xtal.atoms.n) - for a = 1:xtal.atoms.n + for a in 1:(xtal.atoms.n) q = 0.0 if xtal.charges.n != 0 q = xtal.charges.q[a] end - @printf(f, "\n%d %s %f %f %f 0 0 0 0 0 0 0 0 %f", a, xtal.atoms.species[a], xtal.atoms.coords.xf[:, a]..., q, ) + @printf( + f, + "\n%d %s %f %f %f 0 0 0 0 0 0 0 0 %f", + a, + xtal.atoms.species[a], + xtal.atoms.coords.xf[:, a]..., + q, + ) end end if !quiet @@ -930,8 +1076,9 @@ function write_cssr(xtal::Crystal, filename::String; quiet::Bool=false) end end -write_cssr(crystal::Crystal; kwargs...) = write_cssr(crystal, String(split(crystal.name, ".")[1])) - +function write_cssr(crystal::Crystal; kwargs...) + return write_cssr(crystal, String(split(crystal.name, ".")[1])) +end """ crystal = remove_duplicate_atoms_and_charges(crystal, r_tol=0.1, q_tol=0.0001, verbose=false) @@ -939,26 +1086,39 @@ write_cssr(crystal::Crystal; kwargs...) = write_cssr(crystal, String(split(cryst Remove duplicate atoms and charges from a crystal. See [`remove_duplicates`](@ref). # arguments -- `crystal::Crystal`: a crystal -- `r_tol::Float64`: atoms/charges are overlapping if within `r_tol` distance (PBC applied) -- `q_tol::Float64`: charges have the same charge value if their charges are within `q_tol` of each other -- `verbose::Bool`: print off which dupicates are removed + + - `crystal::Crystal`: a crystal + - `r_tol::Float64`: atoms/charges are overlapping if within `r_tol` distance (PBC applied) + - `q_tol::Float64`: charges have the same charge value if their charges are within `q_tol` of each other + - `verbose::Bool`: print off which dupicates are removed # returns -- `crystal::Crystal`: crystal with duplicate atoms and charges removed. + + - `crystal::Crystal`: crystal with duplicate atoms and charges removed. """ -function remove_duplicate_atoms_and_charges(crystal::Crystal, r_tol::Float64=0.1, q_tol::Float64=0.0001) +function remove_duplicate_atoms_and_charges( + crystal::Crystal, + r_tol::Float64=0.1, + q_tol::Float64=0.0001 +) if ne(crystal.bonds) != 0 error("cannot remove duplicates with bonds assigned") end - atoms = remove_duplicates(crystal.atoms, crystal.box, true, r_tol=r_tol, q_tol=q_tol) - charges = remove_duplicates(crystal.charges, crystal.box, true, r_tol=r_tol, q_tol=q_tol) - return Crystal(crystal.name, crystal.box, atoms, charges, MetaGraph(atoms.n), crystal.symmetry) + atoms = remove_duplicates(crystal.atoms, crystal.box, true; r_tol=r_tol, q_tol=q_tol) + charges = + remove_duplicates(crystal.charges, crystal.box, true; r_tol=r_tol, q_tol=q_tol) + return Crystal( + crystal.name, + crystal.box, + atoms, + charges, + MetaGraph(atoms.n), + crystal.symmetry + ) end inside(crystal::Crystal) = inside(crystal.atoms.coords) && inside(crystal.charges.coords) - """ crystal_with_charges = assign_charges(crystal, species_to_charge, net_charge_tol=1e-5) @@ -980,33 +1140,53 @@ crystal_with_charges = assign_charges(crystal, species_to_charge) # tol 1e-5 def ``` # Arguments -- `crystal::Crystal`: the crystal -- `species_to_charge::Dict{Symbol, Float64}`: a dictionary that maps atomic species to charge -- `net_charge_tol::Float64`: the net charge tolerated when asserting charge neutrality of -the resulting crystal + + - `crystal::Crystal`: the crystal + - `species_to_charge::Dict{Symbol, Float64}`: a dictionary that maps atomic species to charge + - `net_charge_tol::Float64`: the net charge tolerated when asserting charge neutrality of + the resulting crystal """ -function assign_charges(crystal::Crystal, species_to_charge::Dict{Symbol, Float64}, net_charge_tol::Float64=1e-5) +function assign_charges( + crystal::Crystal, + species_to_charge::Dict{Symbol, Float64}, + net_charge_tol::Float64=1e-5 +) if crystal.charges.n != 0 - @warn @sprintf("Charges are already present in %s. Removing the current charges on the - crystal and adding new ones...\n", crystal.name) + @warn @sprintf( + "Charges are already present in %s. Removing the current charges on the +crystal and adding new ones...\n", + crystal.name + ) end q = [species_to_charge[s] for s in crystal.atoms.species] charges = Charges(q, crystal.atoms.coords) - new_crystal = Crystal(crystal.name, crystal.box, crystal.atoms, charges, crystal.bonds, crystal.symmetry) + new_crystal = Crystal( + crystal.name, + crystal.box, + crystal.atoms, + charges, + crystal.bonds, + crystal.symmetry + ) # check for charge neutrality - if ! neutral(new_crystal, net_charge_tol) - error(@sprintf("Net charge of crystal %s = %f > net charge tolerance %f. If - charge neutrality is not a problem, pass `net_charge_tol=Inf`\n", crystal.name, - net_charge(new_crystal), net_charge_tol)) + if !neutral(new_crystal, net_charge_tol) + error( + @sprintf( + "Net charge of crystal %s = %f > net charge tolerance %f. If + charge neutrality is not a problem, pass `net_charge_tol=Inf`\n", + crystal.name, + net_charge(new_crystal), + net_charge_tol + ) + ) end return new_crystal end - function Base.show(io::IO, crystal::Crystal) println(io, "Name: ", crystal.name) println(io, crystal.box) @@ -1020,29 +1200,30 @@ function Base.show(io::IO, crystal::Crystal) end @printf(io, "\tbonding graph:\n") println(io, "\t\t# vertices = ", nv(crystal.bonds)) - println(io, "\t\t# edges = ", ne(crystal.bonds)) + return println(io, "\t\t# edges = ", ne(crystal.bonds)) end - function Base.isapprox(c1::Crystal, c2::Crystal; atol::Real=0.0) # name not included - box_flag = isapprox(c1.box, c2.box, atol=atol) + box_flag = isapprox(c1.box, c2.box; atol=atol) if c1.charges.n != c2.charges.n return false end if c1.atoms.n != c2.atoms.n return false end - charges_flag = isapprox(c1.charges, c2.charges, atol=atol) - atoms_flag = isapprox(c1.atoms, c2.atoms, atol=atol) + charges_flag = isapprox(c1.charges, c2.charges; atol=atol) + atoms_flag = isapprox(c1.atoms, c2.atoms; atol=atol) sym1 = c1.symmetry sym2 = c2.symmetry - symmetry_flag = (sym1.is_p1 == sym2.is_p1) && (sym1.space_group == sym1.space_group) && (sym1.operations == sym2.operations) + symmetry_flag = + (sym1.is_p1 == sym2.is_p1) && + (sym1.space_group == sym1.space_group) && + (sym1.operations == sym2.operations) bonds_flag = c1.bonds == c2.bonds return box_flag && charges_flag && atoms_flag && symmetry_flag && bonds_flag end - # slicing of crystal by arrays of Int's function Base.getindex(crystal::Crystal, ids::Union{Array{Int, 1}, UnitRange{Int}}) # mapping from old index to new index @@ -1053,46 +1234,69 @@ function Base.getindex(crystal::Crystal, ids::Union{Array{Int, 1}, UnitRange{Int for edge in edges(crystal.bonds) if edge.src in ids && edge.dst in ids add_edge!(bonds, old_to_new[edge.src], old_to_new[edge.dst]) - set_props!(bonds, old_to_new[edge.src], old_to_new[edge.dst], props(crystal.bonds, edge.src, edge.dst)) + set_props!( + bonds, + old_to_new[edge.src], + old_to_new[edge.dst], + props(crystal.bonds, edge.src, edge.dst) + ) end end if crystal.charges.n == 0 - return Crystal(crystal.name, crystal.box, crystal.atoms[ids], crystal.charges, bonds, crystal.symmetry) + return Crystal( + crystal.name, + crystal.box, + crystal.atoms[ids], + crystal.charges, + bonds, + crystal.symmetry + ) elseif (crystal.charges.n == crystal.atoms.n) && - isapprox(crystal.charges.coords, crystal.atoms.coords) - return Crystal(crystal.name, crystal.box, crystal.atoms[ids], crystal.charges[ids], bonds, crystal.symmetry) + isapprox(crystal.charges.coords, crystal.atoms.coords) + return Crystal( + crystal.name, + crystal.box, + crystal.atoms[ids], + crystal.charges[ids], + bonds, + crystal.symmetry + ) else - error("for getindex(crystal), crystal must have 0 charges or an equal number of charges and atoms that share coordinates") + error( + "for getindex(crystal), crystal must have 0 charges or an equal number of charges and atoms that share coordinates" + ) end return crystal end # slicing of crystal by arrays of Bit's (overload above) -Base.getindex(crystal::Crystal, ids::Union{BitArray{1}}) = getindex(crystal, [i for i = eachindex(ids) if ids[i]]) - +function Base.getindex(crystal::Crystal, ids::Union{BitArray{1}}) + return getindex(crystal, [i for i = eachindex(ids) if ids[i]]) +end function Base.lastindex(crystal::Crystal) if (crystal.atoms.n == crystal.charges.n) || (crystal.charges.n == 0) return crystal.atoms.n else - error("to index the crystal, it must have 0 charges or an equal number of charges and atoms") + error( + "to index the crystal, it must have 0 charges or an equal number of charges and atoms" + ) end end - function Base.:+(crystals::Crystal...; check_overlap::Bool=true, name::String="added_xtal") - box = crystals[1].box + box = crystals[1].box symmetry = crystals[1].symmetry # all crystals must have same boxes and space group - for i = 2:length(crystals) + for i in 2:length(crystals) @assert isapprox(crystals[i].box, box) @assert crystals[i].symmetry.space_group == symmetry.space_group end # initialize atoms, charges, and bonds - atoms = deepcopy(crystals[1].atoms) + atoms = deepcopy(crystals[1].atoms) charges = deepcopy(crystals[1].charges) - bonds = deepcopy(crystals[1].bonds) + bonds = deepcopy(crystals[1].bonds) @assert nv(bonds) == crystals[1].atoms.n for crystal in crystals[2:end] add_vertices!(bonds, crystal.atoms.n) @@ -1100,7 +1304,7 @@ function Base.:+(crystals::Crystal...; check_overlap::Bool=true, name::String="a add_edge!(bonds, atoms.n + ed.src, atoms.n + ed.dst, props(crystal.bonds, ed)) end - atoms += crystal.atoms + atoms += crystal.atoms charges += crystal.charges end @@ -1108,7 +1312,7 @@ function Base.:+(crystals::Crystal...; check_overlap::Bool=true, name::String="a if check_overlap overlap_flag, overlap_pairs = overlap(crystal.atoms.coords, crystal.box, true) - @assert ! overlap_flag "Addition causes overlap: $overlap_pairs" + @assert !overlap_flag "Addition causes overlap: $overlap_pairs" end return crystal @@ -1130,13 +1334,17 @@ atomic_symbol(atom::Atoms) = atom.species[1] function position(crystal::Crystal) pos = Cart(crystal.atoms.coords, crystal.box).x - return [pos[:,i] for i in 1:size(pos,2)] * u"Å" + return [pos[:, i] for i in 1:size(pos, 2)] * u"Å" end velocity(::Crystal) = missing -bounding_box(crystal::Crystal) = SVector{3}([SVector{3}(crystal.box.f_to_c[:,i]u"Å") for i in 1:3]) +function bounding_box(crystal::Crystal) + return SVector{3}([SVector{3}(crystal.box.f_to_c[:, i]u"Å") for i in 1:3]) +end -boundary_conditions(::Crystal) = SVector{3,BoundaryCondition}([Periodic(), Periodic(), Periodic()]) +function boundary_conditions(::Crystal) + return SVector{3, BoundaryCondition}([Periodic(), Periodic(), Periodic()]) +end chemical_formula(crystal::Crystal) = AtomsBase.chemical_formula(crystal) diff --git a/src/distance.jl b/src/distance.jl index bafafec..2f7413c 100644 --- a/src/distance.jl +++ b/src/distance.jl @@ -8,16 +8,15 @@ space; modifies `dxf` for nearest image convention. Fractional coordinates here Warning: this assumes the two molecules are in the box described by fractional coords [0, 1]³. # Arguments -- `dxf::Array{Float64}`: A vector between two atoms in fractional space + + - `dxf::Array{Float64}`: A vector between two atoms in fractional space """ -@inline function nearest_image!(dxf::Array{Float64}) +@inline nearest_image!(dxf::Array{Float64}) = for i in eachindex(dxf) @inbounds if abs(dxf[i]) > 0.5 @inbounds dxf[i] -= sign(dxf[i]) end end -end - @doc raw""" r = distance(coords, box, i, j, apply_pbc) @@ -89,25 +88,29 @@ function distance(coords::Cart, i, j) end distance(atoms::Atoms{Cart}, i, j) = distance(atoms.coords, i, j) -distance(atoms::Atoms, box::Box, i, j, apply_pbc::Bool) = distance(atoms.coords, box, i, j, apply_pbc) -distance(charges::Charges, box::Box, i, j, apply_pbc::Bool) = distance(charges.coords, box, i, j, apply_pbc) - +function distance(atoms::Atoms, box::Box, i, j, apply_pbc::Bool) + return distance(atoms.coords, box, i, j, apply_pbc) +end +function distance(charges::Charges, box::Box, i, j, apply_pbc::Bool) + return distance(charges.coords, box, i, j, apply_pbc) +end function pairwise_distances(coords::Frac, box::Box, apply_pbc::Bool) n = length(coords) pd = zeros(n, n) - for i = 1:n - for j = 1:n + for i in 1:n + for j in 1:n pd[i, j] = distance(coords, box, i, j, apply_pbc) if i > j - pd[i, j] = pd[j, i] + pd[i, j] = pd[j, i] end end end return pd end -pairwise_distances(coords::Cart, box::Box, apply_pbc::Bool) = pairwise_distances(Frac(coords, box), box, apply_pbc) - +function pairwise_distances(coords::Cart, box::Box, apply_pbc::Bool) + return pairwise_distances(Frac(coords, box), box, apply_pbc) +end """ overlap_flag, overlap_pairs = overlap(frac_coords, box, apply_pbc; tol=0.1) @@ -120,23 +123,25 @@ is less than `tol`. Determine if any of the atoms in `crystal` overlap, as determined by comparing their pairwise distances to the lesser of each pair's covalent radii. # Arguments -- `coords::Frac`: the fractional coordinates (`Frac>:Coords`) -- `box::Box`: unit cell information -- `apply_pbc::Bool`: `true` if we wish to apply periodic boundary conditions, `false` otherwise -- `tol::Float64`: tolerance for overlap; if distance between particles less than this, overlap occurs -- `crystal::Crystal`: a crystal to check for overlapping atoms + + - `coords::Frac`: the fractional coordinates (`Frac>:Coords`) + - `box::Box`: unit cell information + - `apply_pbc::Bool`: `true` if we wish to apply periodic boundary conditions, `false` otherwise + - `tol::Float64`: tolerance for overlap; if distance between particles less than this, overlap occurs + - `crystal::Crystal`: a crystal to check for overlapping atoms # Returns -* `overlap_flag::Bool`: `true` if overlap, `false` otherwise -* `overlap_ids::Array{Tuple{Int, Int}, 1}`: ids of coordinate pairs that are overlapping e.g. `[(4, 5), (7, 8)]` + + - `overlap_flag::Bool`: `true` if overlap, `false` otherwise + - `overlap_ids::Array{Tuple{Int, Int}, 1}`: ids of coordinate pairs that are overlapping e.g. `[(4, 5), (7, 8)]` """ function overlap(coords::Frac, box::Box, apply_pbc::Bool; tol::Float64=0.1) overlap_flag = false overlap_ids = Array{Tuple{Int, Int}, 1}() n = size(coords.xf)[2] # number of coords - for i = 1:n - for j = i+1:n + for i in 1:n + for j in (i + 1):n r = distance(coords, box, i, j, apply_pbc) if r < tol push!(overlap_ids, (i, j)) @@ -152,8 +157,8 @@ function overlap(xtal::Crystal, apply_pbc::Bool) overlap_ids = Array{Tuple{Int, Int}, 1}() n = xtal.atoms.n - for i = 1:n - for j = i+1:n + for i in 1:n + for j in (i + 1):n r = distance(xtal.atoms.coords, xtal.box, i, j, apply_pbc) atom_i = strip_number_from_label(xtal.atoms.species[i]) atom_j = strip_number_from_label(xtal.atoms.species[j]) @@ -168,7 +173,6 @@ function overlap(xtal::Crystal, apply_pbc::Bool) return overlap_flag, overlap_ids end - """ atoms = remove_duplicates(atoms, box, apply_pbc, r_tol=0.1) charges = remove_duplicates(charges, box, apply_pbc, r_tol=0.1) @@ -178,25 +182,32 @@ remove duplicates from atoms or charges. loops through all pairs of atoms/charges. if a pair is a duplicate, one is deleted. two atoms are duplicates if both: -* same species -* less than a distance `r_tol` apart -two charges are duplicate if both: -* charge values are within `q_tol` -* less than a distance `r_tol` apart + + - same species + - less than a distance `r_tol` apart + two charges are duplicate if both: + - charge values are within `q_tol` + - less than a distance `r_tol` apart # arguments -- `atoms::Atoms`: the atoms -- `charges::Charges`: the charges -- `box::Box`: unit cell information -- `apply_pbc::Bool`: true iff we apply periodic boundary conditions when computing the distance. -- `r_tol::Float64`: atoms/charges are overlapping if within `r_tol` distance (PBC applied) -- `q_tol::Float64`: charges have the same charge value if their charges are within `q_tol` of each other + + - `atoms::Atoms`: the atoms + - `charges::Charges`: the charges + - `box::Box`: unit cell information + - `apply_pbc::Bool`: true iff we apply periodic boundary conditions when computing the distance. + - `r_tol::Float64`: atoms/charges are overlapping if within `r_tol` distance (PBC applied) + - `q_tol::Float64`: charges have the same charge value if their charges are within `q_tol` of each other """ -function remove_duplicates(ac::Union{Atoms{Frac}, Charges{Frac}}, box::Box, apply_pbc::Bool; - r_tol::Float64=0.1, q_tol::Float64=0.0001) +function remove_duplicates( + ac::Union{Atoms{Frac}, Charges{Frac}}, + box::Box, + apply_pbc::Bool; + r_tol::Float64=0.1, + q_tol::Float64=0.0001 +) ids_keep = trues(ac.n) - for i = 1:ac.n - for j = i+1:ac.n + for i in 1:(ac.n) + for j in (i + 1):(ac.n) if isa(ac, Atoms{Frac}) # if different species, not duplicates. if ac.species[i] != ac.species[j] @@ -204,7 +215,7 @@ function remove_duplicates(ac::Union{Atoms{Frac}, Charges{Frac}}, box::Box, appl end elseif isa(ac, Charges{Frac}) # if different value of charge, not duplicates. - if ! isapprox(ac.q[i], ac.q[j], atol=q_tol) + if !isapprox(ac.q[i], ac.q[j]; atol=q_tol) continue end end @@ -218,7 +229,6 @@ function remove_duplicates(ac::Union{Atoms{Frac}, Charges{Frac}}, box::Box, appl return ac[ids_keep] end - """ dm = distance_matrix(crystal, apply_pbc) @@ -227,16 +237,18 @@ distance between atom `i` and `j`. This matrix is symmetric. If `apply_pbc = tru periodic boundary conditions are applied when computing the distance. # Arguments + -`crystal::Crystal`: crystal structure -`apply_pbc::Bool`: whether or not to apply periodic boundary conditions when computing the distance # Returns + -`dm::Array{Float64, 2}`: symmetric, square distance matrix with zeros on the diagonal """ function distance_matrix(crystal::Crystal, apply_pbc::Bool) dm = zeros(crystal.atoms.n, crystal.atoms.n) - for i = 1:crystal.atoms.n - for j = (i+1):crystal.atoms.n + for i in 1:(crystal.atoms.n) + for j in (i + 1):(crystal.atoms.n) dm[i, j] = distance(crystal.atoms, crystal.box, i, j, apply_pbc) dm[j, i] = dm[i, j] # symmetric end diff --git a/src/matter.jl b/src/matter.jl index c37d82a..b2475af 100644 --- a/src/matter.jl +++ b/src/matter.jl @@ -15,15 +15,16 @@ generally, fractional coordinates should be in [0, 1] and are implicitly associated with a `Box` to represent a periodic coordinate system. e.g. + ```julia f_coords = Frac(rand(3, 2)) # 2 particles f_coords.xf # retreive fractional coords ``` """ -struct Frac<:Coords +struct Frac <: Coords xf::Array{Float64, 2} end -Base.isapprox(c1::Frac, c2::Frac; atol::Real=0) = isapprox(c1.xf, c2.xf, atol=atol) +Base.isapprox(c1::Frac, c2::Frac; atol::Real=0) = isapprox(c1.xf, c2.xf; atol=atol) Base.hcat(c1::Frac, c2::Frac) = Frac(hcat(c1.xf, c2.xf)) """ @@ -32,15 +33,16 @@ cartesian coordinates, a subtype of `Coords`. construct by passing an `Array{Float64, 2}` whose columns are the coordinates. e.g. + ```julia c_coords = Cart(rand(3, 2)) # 2 particles c_coords.x # retreive cartesian coords ``` """ -struct Cart<:Coords +struct Cart <: Coords x::Array{Float64, 2} end -Base.isapprox(c1::Cart, c2::Cart; atol::Real=0) = isapprox(c1.x, c2.x, atol=atol) +Base.isapprox(c1::Cart, c2::Cart; atol::Real=0) = isapprox(c1.x, c2.x; atol=atol) Base.hcat(c1::Cart, c2::Cart) = Cart(hcat(c1.x, c2.x)) # helpers for single vectors @@ -60,7 +62,6 @@ Base.setindex!(coords::Cart, val::Array{Float64, 1}, ids) = (coords.x[:, ids] = Base.size(coords::Cart) = size(coords.x) Base.size(coords::Frac) = size(coords.xf) - """ translate_by!(coords, dx) translate_by!(coords, dx, box) @@ -80,14 +81,9 @@ note that periodic boundary conditions are *not* subsequently applied here. if applied to a `molecule::Molecule`, the coords of atoms, charges, and center of mass are all translated. """ -function translate_by!(coords::Cart, dx::Cart) - coords.x .= broadcast(+, coords.x, dx.x) -end - -function translate_by!(coords::Frac, dxf::Frac) - coords.xf .= broadcast(+, coords.xf, dxf.xf) -end +translate_by!(coords::Cart, dx::Cart) = coords.x .= broadcast(+, coords.x, dx.x) +translate_by!(coords::Frac, dxf::Frac) = coords.xf .= broadcast(+, coords.xf, dxf.xf) """ wrap!(f::Frac) @@ -96,10 +92,7 @@ end wrap fractional coordinates to [0, 1] via `mod(⋅, 1.0)`. e.g. -0.1 --> 0.9 and 1.1 -> 0.1 """ -function wrap!(f::Frac) - f.xf .= mod.(f.xf, 1.0) -end - +wrap!(f::Frac) = f.xf .= mod.(f.xf, 1.0) ### # Atoms @@ -109,7 +102,7 @@ end used to represent a set of atoms in space (their atomic species and coordinates). ```julia -struct Atoms{T<:Coords} # enforce that the type specified is `Coords` +struct Atoms{T <: Coords} # enforce that the type specified is `Coords` n::Int # how many atoms? species::Array{Symbol, 1} # list of species coords::T # coordinates @@ -119,13 +112,14 @@ end here, `T` is `Frac` or `Cart`. helper constructor (infers `n`): + ```julia species = [:H, :H] coords = Cart(rand(3, 2)) atoms = Atoms(species, coords) ``` """ -struct Atoms{T<:Coords} # enforce that the type specified is `Coords` +struct Atoms{T <: Coords} # enforce that the type specified is `Coords` n::Int # how many atoms? species::Array{Symbol, 1} # list of species coords::T # coordinates @@ -138,16 +132,19 @@ function Atoms(species::Array{Symbol, 1}, coords::Coords) end Atoms(species::Symbol, coords::Coords) = Atoms([species], coords) -Atoms{Frac}(n::Int) = Atoms([:_ for a = 1:n], Frac([NaN for i = 1:3, a = 1:n])) # safe pre-allocation -Atoms{Cart}(n::Int) = Atoms([:_ for a = 1:n], Cart([NaN for i = 1:3, a = 1:n])) # safe pre-allocation +Atoms{Frac}(n::Int) = Atoms([:_ for a in 1:n], Frac([NaN for i in 1:3, a in 1:n])) # safe pre-allocation +Atoms{Cart}(n::Int) = Atoms([:_ for a in 1:n], Cart([NaN for i in 1:3, a in 1:n])) # safe pre-allocation -Base.isapprox(a1::Atoms, a2::Atoms; atol::Real=0) = (a1.species == a2.species) && isapprox(a1.coords, a2.coords, atol=atol) -Base.:+(a1::Atoms, a2::Atoms) = Atoms(a1.n + a2.n, [a1.species; a2.species], hcat(a1.coords, a2.coords)) +function Base.isapprox(a1::Atoms, a2::Atoms; atol::Real=0) + return (a1.species == a2.species) && isapprox(a1.coords, a2.coords; atol=atol) +end +function Base.:+(a1::Atoms, a2::Atoms) + return Atoms(a1.n + a2.n, [a1.species; a2.species], hcat(a1.coords, a2.coords)) +end Base.getindex(atoms::Atoms, ids) = Atoms(atoms.species[ids], atoms.coords[ids]) Base.lastindex(atoms::Atoms) = atoms.n - ### # point charges ### @@ -155,7 +152,7 @@ Base.lastindex(atoms::Atoms) = atoms.n used to represent a set of partial point charges in space (their charges and coordinates). ```julia -struct Charges{T<:Coords} # enforce that the type specified is `Coords` +struct Charges{T <: Coords} # enforce that the type specified is `Coords` n::Int q::Array{Float64, 1} coords::T @@ -165,13 +162,14 @@ end here, `T` is `Frac` or `Cart`. helper constructor (infers `n`): + ```julia q = [0.1, -0.1] coords = Cart(rand(3, 2)) charges = Charges(q, coords) ``` """ -struct Charges{T<:Coords} # enforce that the type specified is `Coords` +struct Charges{T <: Coords} # enforce that the type specified is `Coords` n::Int q::Array{Float64, 1} coords::T @@ -183,16 +181,19 @@ function Charges(q::Array{Float64, 1}, coords::Coords) end Charges(q::Float64, coords::Coords) = Charges([q], coords) # for one charge -Charges{Frac}(n::Int) = Charges([NaN for c = 1:n], Frac([NaN for i = 1:3, c = 1:n])) # safe pre-allocation -Charges{Cart}(n::Int) = Charges([NaN for c = 1:n], Cart([NaN for i = 1:3, c = 1:n])) # safe pre-allocation +Charges{Frac}(n::Int) = Charges([NaN for c in 1:n], Frac([NaN for i in 1:3, c in 1:n])) # safe pre-allocation +Charges{Cart}(n::Int) = Charges([NaN for c in 1:n], Cart([NaN for i in 1:3, c in 1:n])) # safe pre-allocation -Base.isapprox(c1::Charges, c2::Charges; atol::Real=0) = isapprox(c1.q, c2.q) && isapprox(c1.coords, c2.coords, atol=atol) -Base.:+(c1::Charges, c2::Charges) = Charges(c1.n + c2.n, [c1.q; c2.q], hcat(c1.coords, c2.coords)) +function Base.isapprox(c1::Charges, c2::Charges; atol::Real=0) + return isapprox(c1.q, c2.q) && isapprox(c1.coords, c2.coords; atol=atol) +end +function Base.:+(c1::Charges, c2::Charges) + return Charges(c1.n + c2.n, [c1.q; c2.q], hcat(c1.coords, c2.coords)) +end Base.getindex(charges::Charges, ids) = Charges(charges.q[ids], charges.coords[ids]) Base.lastindex(charges::Charges) = charges.n - """ nc = net_charge(charges) nc = net_charge(crystal) @@ -201,14 +202,12 @@ Base.lastindex(charges::Charges) = charges.n find the sum of charges in `charges::Charges` or charges in `crystal::Crystal` or `molecule::Molecule`. (if there are no charges, the net charge is zero.) """ -function net_charge(charges::Charges) +net_charge(charges::Charges) = if charges.n == 0 return 0.0 else return sum(charges.q) end -end - """ neutral(charges, tol) # true or false. default tol = 1e-5 diff --git a/src/misc.jl b/src/misc.jl index 8f7f60f..9cb1d1d 100644 --- a/src/misc.jl +++ b/src/misc.jl @@ -1,26 +1,24 @@ -function banner() - FIGlet.render("Xtals.jl", FIGlet.availablefonts()[5]) -end - +banner() = FIGlet.render("Xtals.jl", FIGlet.availablefonts()[5]) function add_extension(filename::String, extension::String) - if ! occursin(extension, filename) + if !occursin(extension, filename) filename *= extension end return filename end - """ atoms = read_xyz("molecule.xyz") read a list of atomic species and their corresponding coordinates from an .xyz file. # Arguments -- `filename::AbstractString`: the path to and filename of the .xyz file + + - `filename::AbstractString`: the path to and filename of the .xyz file # Returns -- `atoms::Atoms{Cart}`: the set of atoms read from the .xyz file. + + - `atoms::Atoms{Cart}`: the set of atoms read from the .xyz file. """ function read_xyz(filename::AbstractString) lines = readlines(filename) @@ -28,16 +26,15 @@ function read_xyz(filename::AbstractString) n = parse(Int, lines[1]) # get number of atoms species = Symbol[] x = zeros(Float64, 3, n) - for i = 1:n + for i in 1:n push!(species, Symbol(split(lines[i + 2])[1])) - for j = 1:3 + for j in 1:3 x[j, i] = parse(Float64, split(lines[i + 2])[1 + j]) end end return Atoms(species, Cart(x)) end - """ write_xyz(atoms, filename; comment="") write_xyz(crystal; comment="", center_at_origin=false) @@ -47,26 +44,36 @@ end write atoms to an .xyz file. # Arguments -- `atoms::Atoms`: the set of atoms. -- `filename::AbstractString`: the filename (absolute path) of the .xyz file. (".xyz" appended automatically -if the extension is not provided.) -- `comment::AbstractString`: comment if you'd like to write to the file. -- `center_at_origin::Bool`: (for crystal only) if `true`, translate all coords such that the origin is the center of the unit cell. + + - `atoms::Atoms`: the set of atoms. + - `filename::AbstractString`: the filename (absolute path) of the .xyz file. (".xyz" appended automatically + if the extension is not provided.) + - `comment::AbstractString`: comment if you'd like to write to the file. + - `center_at_origin::Bool`: (for crystal only) if `true`, translate all coords such that the origin is the center of the unit cell. """ -function write_xyz(atoms::Atoms{Cart}, filename::AbstractString; comment::AbstractString="# $filename") +function write_xyz( + atoms::Atoms{Cart}, + filename::AbstractString; + comment::AbstractString="# $filename" +) filename = add_extension(filename, ".xyz") open(filename, "w") do xyzfile @printf(xyzfile, "%d\n#%s\n", atoms.n, comment) - for i = 1:atoms.n - @printf(xyzfile, "%s %.4f %.4f %.4f\n", atoms.species[i], - atoms.coords.x[1, i], atoms.coords.x[2, i], atoms.coords.x[3, i]) + for i in 1:(atoms.n) + @printf( + xyzfile, + "%s %.4f %.4f %.4f\n", + atoms.species[i], + atoms.coords.x[1, i], + atoms.coords.x[2, i], + atoms.coords.x[3, i] + ) end end return nothing end - """ atoms, bonds = read_mol("molecule.mol") @@ -74,12 +81,14 @@ read a `.mol` file, which contains info about both atoms and bonds. see [here](https://chem.libretexts.org/Courses/University_of_Arkansas_Little_Rock/ChemInformatics_(2017)%3A_Chem_4399%2F%2F5399/2.2%3A_Chemical_Representations_on_Computer%3A_Part_II/2.2.2%3A_Anatomy_of_a_MOL_file) for the anatomy of a `.mol` file. # Arguments -- `filename::AbstractString`: the path to and filename of the .mol file (must pass extension) + + - `filename::AbstractString`: the path to and filename of the .mol file (must pass extension) # Returns -- `atoms::Atoms{Cart}`: the set of atoms read from the `.mol` file. -- `bonds::MetaGraph`: the bonding graph of the atoms read from the `.mol` file. -- `bond_types::Array{Int, 1}`: the array of bond types. + + - `atoms::Atoms{Cart}`: the set of atoms read from the `.mol` file. + - `bonds::MetaGraph`: the bonding graph of the atoms read from the `.mol` file. + - `bond_types::Array{Int, 1}`: the array of bond types. """ function read_mol(filename::String) lines = readlines(filename) @@ -87,21 +96,19 @@ function read_mol(filename::String) n_atoms = parse(Int, split(lines[4])[1]) n_bonds = parse(Int, split(lines[4])[2]) - atoms = Atoms([:blah for i = 1:n_atoms], - Cart(zeros(Float64, 3, n_atoms)) - ) - for i = 1:n_atoms - line_atom_i = split(lines[4+i]) + atoms = Atoms([:blah for i in 1:n_atoms], Cart(zeros(Float64, 3, n_atoms))) + for i in 1:n_atoms + line_atom_i = split(lines[4 + i]) atoms.species[i] = Symbol(line_atom_i[4]) - for j = 1:3 + for j in 1:3 atoms.coords.x[j, i] = parse(Float64, line_atom_i[j]) end end bonds = MetaGraph(n_atoms) - bond_types = [-1 for _ ∈ 1:n_bonds] - for b = 1:n_bonds + bond_types = [-1 for _ in 1:n_bonds] + for b in 1:n_bonds line_bond_b = split(lines[n_atoms + 4 + b]) i = parse(Int, line_bond_b[1]) @@ -114,15 +121,15 @@ function read_mol(filename::String) return atoms, bonds, bond_types end - """ write_mol2(xtal, filename="my_xtal.mol2") Write a `Crystal` to disk in the mol2 format. Includes atoms, bonds, and unit cell. # Arguments -- `xtal::Crystal` : the crystal to export -- `filename::String` : (Optional) the name of the file to save to. By default, file is named automatically from `xtal.name` + + - `xtal::Crystal` : the crystal to export + - `filename::String` : (Optional) the name of the file to save to. By default, file is named automatically from `xtal.name` """ function write_mol2(xtal::Crystal; filename::String=split(xtal.name, ".cif")[1] * ".mol2") # open buffer @@ -134,7 +141,7 @@ function write_mol2(xtal::Crystal; filename::String=split(xtal.name, ".cif")[1] # now the ATOM record @printf(f, "@ATOM\n") coords = Cart(xtal.atoms.coords, xtal.box) - for i in 1:xtal.atoms.n + for i in 1:(xtal.atoms.n) @printf(f, "%d X %f %f %f %s\n", i, coords.x[:, i]..., xtal.atoms.species[i]) end # BOND record @@ -147,25 +154,34 @@ function write_mol2(xtal::Crystal; filename::String=split(xtal.name, ".cif")[1] # unit cell (CRYSIN) @assert xtal.symmetry.is_p1 "Crystal must be in P1 symmetry." @printf(f, "\n@CRYSIN\n") - @printf(f, "%f %f %f %f %f %f 1 1\n", xtal.box.a, xtal.box.b, xtal.box.c, - xtal.box.α, xtal.box.β, xtal.box.γ) + @printf( + f, + "%f %f %f %f %f %f 1 1\n", + xtal.box.a, + xtal.box.b, + xtal.box.c, + xtal.box.α, + xtal.box.β, + xtal.box.γ + ) end end - """ `set_paths("path_to_data", print_paths=true)` Sets all paths in `rc[:paths]` relative to `path_to_data`. Paths follow the standard format of `rc[:paths][:foo] = "path_to_data/foo"`, except for `rc[:paths][:data]` which is `"path_to_data"`. Warnings are issued if any chosen paths are not valid folders. + # Arguments -- `path_to_data::String` : an absolute or relative path to use as the root of the data folder tree. Defaults to present working directory. -- `print_paths::Bool` : Optional. If `true`, prints contents of `rc[:paths]` to console. Defaults to `false`. -- `no_warn::Bool` : Optional. Set `true` to suppress invalid path warnings. Default to `false`. + + - `path_to_data::String` : an absolute or relative path to use as the root of the data folder tree. Defaults to present working directory. + - `print_paths::Bool` : Optional. If `true`, prints contents of `rc[:paths]` to console. Defaults to `false`. + - `no_warn::Bool` : Optional. Set `true` to suppress invalid path warnings. Default to `false`. """ function set_paths(path_to_data::String=pwd(); print_paths::Bool=false, no_warn::Bool=false) - for (key, path) ∈ rc[:paths] # set all relative paths + for (key, path) in rc[:paths] # set all relative paths rc[:paths][key] = joinpath(path_to_data, String(key)) end rc[:paths][:data] = path_to_data # correct data root path @@ -173,7 +189,7 @@ function set_paths(path_to_data::String=pwd(); print_paths::Bool=false, no_warn: @info rc[:paths] end if !no_warn - for (key, path) ∈ rc[:paths] + for (key, path) in rc[:paths] if !isdir(path) @warn "$key path directory not found" path end @@ -181,13 +197,15 @@ function set_paths(path_to_data::String=pwd(); print_paths::Bool=false, no_warn: end end - """ `view_crystal(xtal, drop_cross_pb_bonds=true)` + Launch a GUI window displaying the crystal. + # Arguments -- `xtal::Crystal` : the crystal to display -- `drop_cross_pb_bonds::Bool` : Optional. Set to `true` to remove bonds that extend across periodic boundaries of the unit cell (these tend to mess up the visualization). Defaults to `true`. + + - `xtal::Crystal` : the crystal to display + - `drop_cross_pb_bonds::Bool` : Optional. Set to `true` to remove bonds that extend across periodic boundaries of the unit cell (these tend to mess up the visualization). Defaults to `true`. """ function view_crystal(xtal::Crystal; drop_cross_pb_bonds::Bool=true) structure_file = tempname() * ".mol2" @@ -196,9 +214,9 @@ function view_crystal(xtal::Crystal; drop_cross_pb_bonds::Bool=true) if drop_cross_pb_bonds drop_cross_pb_bonds!(crystal) end - write_mol2(crystal, filename=structure_file) + write_mol2(crystal; filename=structure_file) write_vtk(crystal.box, box_file) - viewfile(structure_file, "mol2", vtkcell=box_file) + viewfile(structure_file, "mol2"; vtkcell=box_file) rm(structure_file) - rm(box_file) + return rm(box_file) end diff --git a/src/repfactors.jl b/src/repfactors.jl index 6f690cf..cd8f16a 100644 --- a/src/repfactors.jl +++ b/src/repfactors.jl @@ -1,15 +1,20 @@ """ repfactors = replication_factors(unitcell, cutoffradius) + Find the replication factors needed to make a supercell big enough to fit a sphere with the specified cutoff radius. In Xtals.jl, rather than replicating the atoms in the home unit cell to build the supercell that serves as a simulation box, we replicate the home unit cell to form the supercell (simulation box) in a for loop. This function ensures enough replication factors such that the nearest image convention can be applied. A non-replicated supercell has 1 as the replication factor in each dimension (`repfactors = (1, 1, 1)`). + # Arguments -- `unitcell::Box`: The unit cell of the framework -- `cutoff_radius::Float64`: Cutoff radius beyond which we define the potential energy to be zero (units: Angstrom) -# Returns -- `repfactors::Tuple{Int, Int, Int}`: The replication factors in the a, b, c directions + + - `unitcell::Box`: The unit cell of the framework + - `cutoff_radius::Float64`: Cutoff radius beyond which we define the potential energy to be zero (units: Angstrom) + +# Returns # Unit vectors used to transform from fractional coordinates to cartesian coordinates. We'll be + + - `repfactors::Tuple{Int, Int, Int}`: The replication factors in the a, b, c directions """ function replication_factors(unitcell::Box, cutoff_radius::Float64) # Unit vectors used to transform from fractional coordinates to cartesian coordinates. We'll be @@ -22,7 +27,7 @@ function replication_factors(unitcell::Box, cutoff_radius::Float64) n_bc = cross(b, c) # c0 defines a center in the unit cell - c0 = [a b c] * [.5, .5, .5] + c0 = [a b c] * [0.5, 0.5, 0.5] rep = [1, 1, 1] @@ -30,22 +35,22 @@ function replication_factors(unitcell::Box, cutoff_radius::Float64) # |n_bc ⋅ c0|/|n_bc| defines the distance from the end of the supercell and the center. As long as that distance is less than the cutoff radius, we need to increase it while abs(dot(n_bc, c0)) / norm(n_bc) < cutoff_radius rep[1] += 1 - a += unitcell.f_to_c[:,1] - c0 = [a b c] * [.5, .5, .5] + a += unitcell.f_to_c[:, 1] + c0 = [a b c] * [0.5, 0.5, 0.5] end # Repeat for `b` while abs(dot(n_ac, c0)) / norm(n_ac) < cutoff_radius rep[2] += 1 - b += unitcell.f_to_c[:,2] - c0 = [a b c] * [.5, .5, .5] + b += unitcell.f_to_c[:, 2] + c0 = [a b c] * [0.5, 0.5, 0.5] end # Repeat for `c` while abs(dot(n_ab, c0)) / norm(n_ab) < cutoff_radius rep[3] += 1 - c += unitcell.f_to_c[:,3] - c0 = [a b c] * [.5, .5, .5] + c += unitcell.f_to_c[:, 3] + c0 = [a b c] * [0.5, 0.5, 0.5] end return (rep[1], rep[2], rep[3])::Tuple{Int, Int, Int} diff --git a/test/aqua.jl b/test/aqua.jl index 5e1a72c..bf31d60 100644 --- a/test/aqua.jl +++ b/test/aqua.jl @@ -2,13 +2,10 @@ using Xtals import Aqua # ambiguity testing finds many "problems" outside the scope of this package -ambiguities=false +ambiguities = false # to skip when checking for stale dependencies and missing compat entries # Aqua is added in a separate CI job, so (ironically) does not work w/ itself stale_deps = (ignore=[:Aqua],) -Aqua.test_all(Xtals; - ambiguities=ambiguities, - stale_deps=stale_deps -) \ No newline at end of file +Aqua.test_all(Xtals; ambiguities=ambiguities, stale_deps=stale_deps) diff --git a/test/assert_p1_symmetry.jl b/test/assert_p1_symmetry.jl index 077979f..7be9ef0 100644 --- a/test/assert_p1_symmetry.jl +++ b/test/assert_p1_symmetry.jl @@ -7,7 +7,7 @@ using IOCapture, Test, Xtals # coefficient test @testset "P1 Symmetry Tests" begin IOCapture.capture() do - non_P1_framework = Crystal("symmetry_test_structure.cif", convert_to_p1=false) + non_P1_framework = Crystal("symmetry_test_structure.cif"; convert_to_p1=false) # Test that an assertion is thrown when trying to replicate a non-P1 # structure @test_throws ErrorException replicate(non_P1_framework, (2, 2, 2)) diff --git a/test/bond_vectors.jl b/test/bond_vectors.jl index 930a54d..9327306 100644 --- a/test/bond_vectors.jl +++ b/test/bond_vectors.jl @@ -3,21 +3,18 @@ xtal = Crystal( "vector test 1", box, - Atoms( - [:H, :H, :H], - Frac([ - 0.000 0.000 0.075; - 0.075 0.000 0.000; + Atoms([:H, :H, :H], Frac([ + 0.000 0.000 0.075 + 0.075 0.000 0.000 0.000 0.000 0.000 - ]) - ), + ])), Charges{Frac}(0) ) # check non-fatal error for calculating vectors w/o bonds @test isnothing(Xtals.calculate_bond_vectors!(xtal)) - infer_bonds!(xtal, true, calculate_vectors=true) + infer_bonds!(xtal, true; calculate_vectors=true) # check fatal error for attempting to re-calculate vectors @test_throws ErrorException calculate_bond_vectors!(xtal) @@ -32,12 +29,12 @@ @test get_prop(xtal.bonds, :has_vectors) # check vectors - @test isapprox(get_bond_vector(xtal, 1, 2), [0.00; -0.75; 0.00]) + @test isapprox(get_bond_vector(xtal, 1, 2), [0.00; -0.75; 0.00]) - @test isapprox(get_bond_vector(xtal, 2, 3), [0.75; 0.00; 0.00]) + @test isapprox(get_bond_vector(xtal, 2, 3), [0.75; 0.00; 0.00]) # check reversal of direction - @test isapprox(get_bond_vector(xtal, 2, 1), [0.00; 0.75; 0.00]) + @test isapprox(get_bond_vector(xtal, 2, 1), [0.00; 0.75; 0.00]) # check bond angle @test isapprox(bond_angle(xtal, 1, 2, 3), deg2rad(90)) @@ -51,79 +48,68 @@ xtal = Crystal( "vector test 2", box, - Atoms( - [:H, :H, :H], - Frac([ - 1.000 0.000 0.075; - 0.075 0.000 0.000; + Atoms([:H, :H, :H], Frac([ + 1.000 0.000 0.075 + 0.075 0.000 0.000 0.000 0.000 0.000 - ]) - ), + ])), Charges{Frac}(0) ) - infer_bonds!(xtal, true, calculate_vectors=true) + infer_bonds!(xtal, true; calculate_vectors=true) + + @test isapprox(get_bond_vector(xtal, 1, 2), [0.00; -0.75; 0.00]) - @test isapprox(get_bond_vector(xtal, 1, 2), [0.00; -0.75; 0.00]) + @test isapprox(get_bond_vector(xtal, 2, 1), [0.00; 0.75; 0.00]) + + @test isapprox(get_bond_vector(xtal, 2, 3), [0.75; 0.00; 0.00]) - @test isapprox(get_bond_vector(xtal, 2, 1), [0.00; 0.75; 0.00]) - - @test isapprox(get_bond_vector(xtal, 2, 3), [0.75; 0.00; 0.00]) - @test isapprox(bond_angle(xtal, 1, 2, 3), deg2rad(90)) # 90° angles work to numerical precision xtal = Crystal( "vector test 3", box, - Atoms( - [:H, :H, :H], - Frac([ - 0.000 0.075 0.150; - 0.000 0.000 0.000; + Atoms([:H, :H, :H], Frac([ + 0.000 0.075 0.150 0.000 0.000 0.000 - ]) - ), + 0.000 0.000 0.000 + ])), Charges{Frac}(0) ) - infer_bonds!(xtal, true, calculate_vectors=true) + infer_bonds!(xtal, true; calculate_vectors=true) @test isapprox(get_bond_vector(xtal, 1, 2), [0.75; 0.00; 0.00]) - + @test isapprox(get_bond_vector(xtal, 2, 3), [0.75; 0.00; 0.00]) - + @test isapprox(bond_angle(xtal, 1, 2, 3), π, rtol=1e-4) # 180° angles work to within 0.01% methane = Frac(read_xyz(joinpath(rc[:paths][:crystals], "methane.xyz")), box) - xtal = Crystal( - "vector test 4", - box, - methane, - Charges{Frac}(0) - ) - infer_bonds!(xtal, true, calculate_vectors=true) - + xtal = Crystal("vector test 4", box, methane, Charges{Frac}(0)) + infer_bonds!(xtal, true; calculate_vectors=true) + # test that ∠ijk = ∠kji and all angles are all the same α = bond_angle(xtal, 2, 1, 3) results = Bool[] - for i ∈ 2:5 - for j ∈ 2:5 + for i in 2:5 + for j in 2:5 if i == j push!(results, isapprox(bond_angle(xtal, i, 1, j), 0)) # 0° angles work to numerical precision else - push!(results, isapprox(bond_angle(xtal, i, 1, j), α, rtol=1e-5)) # identical angles work to within 0.001% + push!(results, isapprox(bond_angle(xtal, i, 1, j), α; rtol=1e-5)) # identical angles work to within 0.001% end - end + end end - + @test all(results) - + # test that bond angles match expected for tetrahedral carbon @test isapprox(bond_angle(xtal, 2, 1, 3), deg2rad(109.5), rtol=1e-3) # 109.5° angles work to within 0.1% - + # test that result is translation invariant results = Vector{Bool}(undef, 10000) - for i ∈ eachindex(results) + for i in eachindex(results) translate_by!(xtal.atoms.coords, Frac(rand(3, 1))) - results[i] = isapprox(bond_angle(xtal, 2, 1, 3), α, rtol=1e-5) + results[i] = isapprox(bond_angle(xtal, 2, 1, 3), α; rtol=1e-5) end - + @test all(results) -end \ No newline at end of file +end diff --git a/test/bonds.jl b/test/bonds.jl index 16083a9..ca95160 100644 --- a/test/bonds.jl +++ b/test/bonds.jl @@ -11,8 +11,8 @@ cof102 = Crystal("cof-102.cif") xtal1 = deepcopy(xtal) infer_bonds!(xtal, true) xtal2 = xtal[1:10] - xtal3 = +(xtal, xtal2, check_overlap=false) - xtal4 = +(xtal2, xtal, check_overlap=false) + xtal3 = +(xtal, xtal2; check_overlap=false) + xtal4 = +(xtal2, xtal; check_overlap=false) # check "conservation of matter" @test nv(xtal3.bonds) == nv(xtal.bonds) + nv(xtal2.bonds) @@ -28,9 +28,9 @@ cof102 = Crystal("cof-102.cif") @test nv(xtal3.bonds) == nv(xtal4.bonds) # check adding when 1 or both have no bonds - xtal5 = +(xtal1, xtal2, check_overlap=false) - xtal6 = +(xtal2, xtal1, check_overlap=false) - xtal7 = +(xtal1, xtal1, check_overlap=false) + xtal5 = +(xtal1, xtal2; check_overlap=false) + xtal6 = +(xtal2, xtal1; check_overlap=false) + xtal7 = +(xtal1, xtal1; check_overlap=false) @test ne(xtal5.bonds) == ne(xtal2.bonds) @@ -44,7 +44,12 @@ include("bond_vectors.jl") @testset "bond sanity check" begin box = unit_cube() texas_carbon = read_xyz(joinpath(rc[:paths][:crystals], "texas_carbon.xyz")) - xtal = Crystal("CH5", box, Atoms(texas_carbon.n, texas_carbon.species, Frac(texas_carbon.coords, box)), Charges{Frac}(0)) + xtal = Crystal( + "CH5", + box, + Atoms(texas_carbon.n, texas_carbon.species, Frac(texas_carbon.coords, box)), + Charges{Frac}(0) + ) @test !bond_sanity_check(xtal)[1] @@ -53,23 +58,29 @@ include("bond_vectors.jl") @test !bond_sanity_check(xtal)[1] trihydrogen = read_xyz(joinpath(rc[:paths][:crystals], "trihydrogen.xyz")) - xtal = Crystal("H3", box, Atoms(trihydrogen.n, trihydrogen.species, Frac(trihydrogen.coords, box)), Charges{Frac}(0)) + xtal = Crystal( + "H3", + box, + Atoms(trihydrogen.n, trihydrogen.species, Frac(trihydrogen.coords, box)), + Charges{Frac}(0) + ) infer_bonds!(xtal, false) @test !bond_sanity_check(xtal)[1] end - @testset "NiPyC2 Tests" begin # NiPyC2 bonding c = Crystal("NiPyC2_relax.cif") strip_numbers_from_atom_labels!(c) c = replicate(c, (4, 4, 4)) - bonding_rules = [BondingRule(:H, :*, 1.2), - BondingRule(:Ni, :O, 2.5), - BondingRule(:Ni, :N, 2.5), - BondingRule(:*, :*, 1.9)] - infer_bonds!(c, true, bonding_rules=bonding_rules) + bonding_rules = [ + BondingRule(:H, :*, 1.2), + BondingRule(:Ni, :O, 2.5), + BondingRule(:Ni, :N, 2.5), + BondingRule(:*, :*, 1.9) + ] + infer_bonds!(c, true; bonding_rules=bonding_rules) c1 = deepcopy(c) conn_comps = connected_components(c.bonds) @@ -83,7 +94,6 @@ end @test length(conn_comps) == 2 # interpenetrated end - @testset "BACMOH Tests" begin xtal = Crystal("BACMOH_clean.cif") strip_numbers_from_atom_labels!(xtal) @@ -96,7 +106,6 @@ end @test xtal.atoms.species[neighbors(xtal.bonds, 1)] == [:O, :N, :N, :O] end - @testset ".mol/.cif bonds vs. inferred" begin mol_atoms, mol_bonds = read_mol("data/example.mol") box = unit_cube() @@ -112,7 +121,6 @@ end @test mol_bonds == xtal.bonds end - @testset "metadata" begin xtal = deepcopy(cof102) infer_bonds!(xtal, true) @@ -129,10 +137,9 @@ end @test xtal.bonds == xtal2.bonds - write_mol2(xtal, filename=tempname()) + write_mol2(xtal; filename=tempname()) end - @testset "bonding rules" begin bonding_rules_path = tempname() write_bonding_rules(bonding_rules_path) @@ -144,12 +151,12 @@ end @test length(bonding_rules) == n - add_bonding_rules([BondingRule(:foo, :bar, 5.)]) + add_bonding_rules([BondingRule(:foo, :bar, 5.0)]) @test length(rc[:bonding_rules]) == (n + 1) capture() do - println([r for r ∈ bonding_rules if :C == r.species_j][1]) + return println([r for r ∈ bonding_rules if :C == r.species_j][1]) end @test true @@ -164,12 +171,11 @@ end @test true end - @testset "etc" begin temp_vtk_path = tempname() xtal = Crystal("SBMOF-1.cif") capture() do - write_bond_information(xtal, temp_vtk_path, center_at_origin=true) + return write_bond_information(xtal, temp_vtk_path; center_at_origin=true) end @test isfile(temp_vtk_path * ".vtk") @@ -181,19 +187,22 @@ end no_pb_temppath = tempname() capture() do write_bond_information(xtal, all_bonds_temppath) - write_bond_information(xtal, no_pb_temppath, bond_filter=:cross_boundary=>p->!p) + return write_bond_information( + xtal, + no_pb_temppath; + bond_filter=:cross_boundary => p -> !p + ) end @test all(isfile.([all_bonds_temppath, no_pb_temppath] .* ".vtk")) capture() do - write_bond_information(xtal; verbose=true) + return write_bond_information(xtal; verbose=true) end @test isfile(xtal.name * "_bonds.vtk") rm(xtal.name * "_bonds.vtk") end - @testset "bonds from xyz" begin xtal = deepcopy(irmof1) infer_bonds!(xtal, false) @@ -202,10 +211,9 @@ end @test bonds == xtal.bonds end - @testset "infer_bonds kwarg" begin - xtal1 = Crystal("IRMOF-1.cif", infer_bonds=true, periodic_boundaries=true) - xtal2 = Crystal("IRMOF-1.cif", infer_bonds=true, periodic_boundaries=false) + xtal1 = Crystal("IRMOF-1.cif"; infer_bonds=true, periodic_boundaries=true) + xtal2 = Crystal("IRMOF-1.cif"; infer_bonds=true, periodic_boundaries=false) @test xtal1.bonds ≠ xtal2.bonds diff --git a/test/box.jl b/test/box.jl index 07bd63a..3751358 100644 --- a/test/box.jl +++ b/test/box.jl @@ -6,8 +6,7 @@ import IOCapture.capture @testset "Box Tests" begin box = Box(11.6, 5.5, 22.9, 90.0 * π / 180, 100.8 * π / 180.0, 90.0 * π / 180.0) @test isapprox(box, Box(box.f_to_c)) # alt. constructor - @test isapprox(box, Box(box.a, box.b, box.c, - box.α, box.β, box.γ)) # alt. constructor + @test isapprox(box, Box(box.a, box.b, box.c, box.α, box.β, box.γ)) # alt. constructor @test box.f_to_c * box.c_to_f ≈ Matrix{Float64}(I, 3, 3) @test isapprox(box.reciprocal_lattice, 2 * π * inv(box.f_to_c)) #TODO just use this method to construct recip. lattice. @test isapprox(replicate(box, (1, 1, 1)), box) @@ -27,10 +26,11 @@ import IOCapture.capture @test isapprox(replicate(box, (3, 5, 4)), Box(3.0, 5.0, 4.0)) box = Box(10.0, 10.0, 10.0) - f = Frac([0.1 0.1; - 0.3 0.2; - 0.5 0.6] - ) + f = Frac([ + 0.1 0.1 + 0.3 0.2 + 0.5 0.6 + ]) atoms = Atoms([:C, :O], f) atoms_c = Cart(atoms, box) @test atoms_c.species == atoms.species @@ -44,20 +44,22 @@ import IOCapture.capture # inside box function box = Box(1.0, 20.0, 30.0) - f = Frac([0.1 0.6; - 0.9 0.8; - 0.8 0.9] - ) + f = Frac([ + 0.1 0.6 + 0.9 0.8 + 0.8 0.9 + ]) @test inside(f) f.xf[2, 2] = -0.1 - @test ! inside(f) + @test !inside(f) f.xf[2, 2] = 1.2 - @test ! inside(f) + @test !inside(f) - c = Cart([0.1 0.4; - 4.5 13.0; - 22. 10.1] - ) + c = Cart([ + 0.1 0.4 + 4.5 13.0 + 22.0 10.1 + ]) @test inside(c, box) c.x[2, 2] = -0.1 @test !inside(c, box) @@ -71,39 +73,43 @@ import IOCapture.capture # translate_by! box = Box(1.0, 10.0, 40.0) - f = Frac([0.1 0.2; - 0.2 0.5; - 0.3 0.8] - ) + f = Frac([ + 0.1 0.2 + 0.2 0.5 + 0.3 0.8 + ]) dx = Cart([0.2, 1.0, 4.0]) translate_by!(f, dx, box) - @test isapprox(f, Frac([0.3 0.4; - 0.3 0.6; - 0.4 0.9] - ) - ) + @test isapprox(f, Frac([ + 0.3 0.4 + 0.3 0.6 + 0.4 0.9 + ])) - c = Cart([1.0 0.0; - 10.0 0.0; - 40.0 0.0]) #corner of box + c = Cart([ + 1.0 0.0 + 10.0 0.0 + 40.0 0.0 + ]) #corner of box dxf = Frac([-1.0, -1.0, -1.0]) translate_by!(c, dxf, box) - @test isapprox(c, Cart([0.0 -1.0; - 0.0 -10.0; - 0.0 -40.0]) - ) + @test isapprox(c, Cart([ + 0.0 -1.0 + 0.0 -10.0 + 0.0 -40.0 + ])) - xtal = Crystal("SBMOF-1.cif") - vtk_temp = tempname() - capture() do - write_vtk(xtal.box, vtk_temp, verbose=true, center_at_origin=true) - end - @test isfile(vtk_temp * ".vtk") + xtal = Crystal("SBMOF-1.cif") + vtk_temp = tempname() + capture() do + return write_vtk(xtal.box, vtk_temp; verbose=true, center_at_origin=true) + end + @test isfile(vtk_temp * ".vtk") - # effective test of Base.show(io::IO, box::Box) - capture() do - println(xtal.box) - end - @test true + # effective test of Base.show(io::IO, box::Box) + capture() do + return println(xtal.box) + end + @test true end end diff --git a/test/crystal.jl b/test/crystal.jl index 4a26018..023e5ec 100644 --- a/test/crystal.jl +++ b/test/crystal.jl @@ -6,9 +6,11 @@ import IOCapture.capture # for test only # if the multi sets are equal, then when you remove duplicates, # you will be left with ac1. -function equal_multisets(ac1::Union{Atoms{Frac}, Charges{Frac}}, - ac2::Union{Atoms{Frac}, Charges{Frac}}, - box::Box) +function equal_multisets( + ac1::Union{Atoms{Frac}, Charges{Frac}}, + ac2::Union{Atoms{Frac}, Charges{Frac}}, + box::Box +) ac = ac1 + ac2 ac_dm = remove_duplicates(ac, box, true) return isapprox(ac_dm, ac1) @@ -16,121 +18,148 @@ end # load here for use in multiple tests sbmof1 = Crystal("SBMOF-1.cif") -sbmof1_no_overlap = Crystal("SBMOF-1_overlap.cif", remove_duplicates=true) +sbmof1_no_overlap = Crystal("SBMOF-1_overlap.cif"; remove_duplicates=true) test_structure2 = Crystal("test_structure2.cif") symmetry_test_structure = Crystal("symmetry_test_structure.cif") symmetry_test_structure_cartn = Crystal("symmetry_test_structure_cartn.cif") irmof1 = Crystal("IRMOF-1.cif") @testset "test data loading" begin - @test !any(isnothing.([sbmof1, sbmof1_no_overlap, test_structure2, symmetry_test_structure, symmetry_test_structure_cartn, irmof1])) + @test !any( + isnothing.([ + sbmof1, + sbmof1_no_overlap, + test_structure2, + symmetry_test_structure, + symmetry_test_structure_cartn, + irmof1 + ]) + ) end @testset "comparisons" begin - # charge/atom count mismatch tests for isapprox + # charge/atom count mismatch tests for isapprox xtal = deepcopy(sbmof1) xtal2 = deepcopy(irmof1) @test !isapprox(xtal, xtal2) # different atoms.n - xtal2 = assign_charges(xtal, Dict(:Ca => -2.0, :O => 2.0, :S => 2.0, :C => 0.0, :H => 0.0), Inf) + xtal2 = assign_charges( + xtal, + Dict(:Ca => -2.0, :O => 2.0, :S => 2.0, :C => 0.0, :H => 0.0), + Inf + ) @test !isapprox(xtal, xtal2) # different charges.n end @testset "symmetry" begin - # test no-op for xtal symmetry when already P1 - xtal1 = deepcopy(sbmof1_no_overlap) + # test no-op for xtal symmetry when already P1 + xtal1 = deepcopy(sbmof1_no_overlap) @test xtal1 == apply_symmetry_operations(xtal1) - # make sure that CIF files with tabs in the _symmetry_equiv_pos_as_xyz blocks load w/o error - xtal = Crystal("xtal_w_tabs.cif", check_overlap=false) + # make sure that CIF files with tabs in the _symmetry_equiv_pos_as_xyz blocks load w/o error + xtal = Crystal("xtal_w_tabs.cif"; check_overlap=false) @test true ##! this test should really be a comparison of xtal vs. something; otherwise, actually only testing the "no error" part - ### apply_symmetry_operations - # test .cif reader for non-P1 symmetry - # no atoms should overlap - # should place atoms in the same positions as the P1 conversion using - # openBabel - # test wraps coords to [0, 1] for symmetry ops - xtal = Crystal("symmetry_test_structure.cif", wrap_coords=false) - @test any(xtal.atoms.coords.xf .> 1.0) - xtal = deepcopy(symmetry_test_structure) # default is to wrap coords - @test all(xtal.atoms.coords.xf .< 1.0) && all(xtal.atoms.coords.xf .> 0.0) - xtal = Crystal("symmetry_test_structure_charges.cif") - @test xtal.charges.n == xtal.atoms.n - @test !any(isnan.(xtal.charges.q)) - - non_P1_crystal = deepcopy(symmetry_test_structure) # not in P1 original - strip_numbers_from_atom_labels!(non_P1_crystal) - - P1_crystal = Crystal("symmetry_test_structure_P1.cif") # in P1 originally - strip_numbers_from_atom_labels!(P1_crystal) - - @test isapprox(non_P1_crystal.box, P1_crystal.box) - @test equal_multisets(non_P1_crystal.atoms, P1_crystal.atoms, P1_crystal.box) - @test equal_multisets(non_P1_crystal.charges, P1_crystal.charges, P1_crystal.box) - - non_P1_cartesian = deepcopy(symmetry_test_structure_cartn) # not in P1 originally - strip_numbers_from_atom_labels!(non_P1_cartesian) - - @test isapprox(non_P1_crystal.box, non_P1_cartesian.box) - @test equal_multisets(non_P1_crystal.atoms, non_P1_cartesian.atoms, P1_crystal.box) - @test equal_multisets(non_P1_crystal.charges, non_P1_cartesian.charges, non_P1_cartesian.box) - - # should all be in P1 now. - @test non_P1_crystal.symmetry.is_p1 - @test non_P1_cartesian.symmetry.is_p1 - @test P1_crystal.symmetry.is_p1 - - # test reading in non-P1 then applying symmetry later - # read in the same files as above, then convert to P1, then compare - non_P1_xtal = Crystal("symmetry_test_structure.cif", convert_to_p1=false) - strip_numbers_from_atom_labels!(non_P1_xtal) - non_P1_xtal_cartesian = Crystal("symmetry_test_structure_cartn.cif", convert_to_p1=false) - strip_numbers_from_atom_labels!(non_P1_xtal_cartesian) - - # make sure these crystals are not in P1 symmetry when convert_to_p1 is - # set to false - @test ! non_P1_xtal.symmetry.is_p1 - @test ! non_P1_xtal_cartesian.symmetry.is_p1 - - # test write_cif in non_p1 symmetry - write_cif(non_P1_xtal, joinpath("data", "crystals", "rewritten_symmetry_test_structure.cif")) - # keep this in cartesian to test both - write_cif(non_P1_xtal_cartesian, joinpath("data", "crystals", "rewritten_symmetry_test_structure_cartn.cif"), fractional_coords=false) - rewritten_non_p1_fractional = Crystal("rewritten_symmetry_test_structure.cif"; convert_to_p1=false) - strip_numbers_from_atom_labels!(rewritten_non_p1_fractional) - rewritten_non_p1_cartesian = Crystal("rewritten_symmetry_test_structure_cartn.cif"; convert_to_p1=false) - strip_numbers_from_atom_labels!(rewritten_non_p1_cartesian) - - @test isapprox(rewritten_non_p1_fractional, non_P1_xtal, atol=0.0001) - @test isapprox(rewritten_non_p1_cartesian, non_P1_xtal_cartesian, atol=0.0001) - - non_P1_xtal = deepcopy(symmetry_test_structure) - strip_numbers_from_atom_labels!(non_P1_xtal) - non_P1_xtal_cartesian = deepcopy(symmetry_test_structure_cartn) - strip_numbers_from_atom_labels!(non_P1_xtal_cartesian) - @test non_P1_xtal.symmetry.is_p1 - @test non_P1_xtal_cartesian.symmetry.is_p1 - - # test that same structure is created when reading and converting to P1 and - # when reading then converting to P1 - @test isapprox(non_P1_crystal, non_P1_xtal) - @test isapprox(non_P1_cartesian, non_P1_xtal_cartesian) + ### apply_symmetry_operations + # test .cif reader for non-P1 symmetry + # no atoms should overlap + # should place atoms in the same positions as the P1 conversion using + # openBabel + # test wraps coords to [0, 1] for symmetry ops + xtal = Crystal("symmetry_test_structure.cif"; wrap_coords=false) + @test any(xtal.atoms.coords.xf .> 1.0) + xtal = deepcopy(symmetry_test_structure) # default is to wrap coords + @test all(xtal.atoms.coords.xf .< 1.0) && all(xtal.atoms.coords.xf .> 0.0) + xtal = Crystal("symmetry_test_structure_charges.cif") + @test xtal.charges.n == xtal.atoms.n + @test !any(isnan.(xtal.charges.q)) + + non_P1_crystal = deepcopy(symmetry_test_structure) # not in P1 original + strip_numbers_from_atom_labels!(non_P1_crystal) + + P1_crystal = Crystal("symmetry_test_structure_P1.cif") # in P1 originally + strip_numbers_from_atom_labels!(P1_crystal) + + @test isapprox(non_P1_crystal.box, P1_crystal.box) + @test equal_multisets(non_P1_crystal.atoms, P1_crystal.atoms, P1_crystal.box) + @test equal_multisets(non_P1_crystal.charges, P1_crystal.charges, P1_crystal.box) + + non_P1_cartesian = deepcopy(symmetry_test_structure_cartn) # not in P1 originally + strip_numbers_from_atom_labels!(non_P1_cartesian) + + @test isapprox(non_P1_crystal.box, non_P1_cartesian.box) + @test equal_multisets(non_P1_crystal.atoms, non_P1_cartesian.atoms, P1_crystal.box) + @test equal_multisets( + non_P1_crystal.charges, + non_P1_cartesian.charges, + non_P1_cartesian.box + ) + + # should all be in P1 now. + @test non_P1_crystal.symmetry.is_p1 + @test non_P1_cartesian.symmetry.is_p1 + @test P1_crystal.symmetry.is_p1 + + # test reading in non-P1 then applying symmetry later + # read in the same files as above, then convert to P1, then compare + non_P1_xtal = Crystal("symmetry_test_structure.cif"; convert_to_p1=false) + strip_numbers_from_atom_labels!(non_P1_xtal) + non_P1_xtal_cartesian = + Crystal("symmetry_test_structure_cartn.cif"; convert_to_p1=false) + strip_numbers_from_atom_labels!(non_P1_xtal_cartesian) + + # make sure these crystals are not in P1 symmetry when convert_to_p1 is + # set to false + @test !non_P1_xtal.symmetry.is_p1 + @test !non_P1_xtal_cartesian.symmetry.is_p1 + + # test write_cif in non_p1 symmetry + write_cif( + non_P1_xtal, + joinpath("data", "crystals", "rewritten_symmetry_test_structure.cif") + ) + # keep this in cartesian to test both + write_cif( + non_P1_xtal_cartesian, + joinpath("data", "crystals", "rewritten_symmetry_test_structure_cartn.cif"); + fractional_coords=false + ) + rewritten_non_p1_fractional = + Crystal("rewritten_symmetry_test_structure.cif"; convert_to_p1=false) + strip_numbers_from_atom_labels!(rewritten_non_p1_fractional) + rewritten_non_p1_cartesian = + Crystal("rewritten_symmetry_test_structure_cartn.cif"; convert_to_p1=false) + strip_numbers_from_atom_labels!(rewritten_non_p1_cartesian) + + @test isapprox(rewritten_non_p1_fractional, non_P1_xtal, atol=0.0001) + @test isapprox(rewritten_non_p1_cartesian, non_P1_xtal_cartesian, atol=0.0001) + + non_P1_xtal = deepcopy(symmetry_test_structure) + strip_numbers_from_atom_labels!(non_P1_xtal) + non_P1_xtal_cartesian = deepcopy(symmetry_test_structure_cartn) + strip_numbers_from_atom_labels!(non_P1_xtal_cartesian) + @test non_P1_xtal.symmetry.is_p1 + @test non_P1_xtal_cartesian.symmetry.is_p1 + + # test that same structure is created when reading and converting to P1 and + # when reading then converting to P1 + @test isapprox(non_P1_crystal, non_P1_xtal) + @test isapprox(non_P1_cartesian, non_P1_xtal_cartesian) end @testset "printing" begin - # test xtal info printing (but don't actually clutter the outputs; print to devnull) + # test xtal info printing (but don't actually clutter the outputs; print to devnull) println(devnull, sbmof1) @test true # test verbose printing in empirical_formula capture() do - empirical_formula(sbmof1, verbose=true) + return empirical_formula(sbmof1; verbose=true) end @test true ##! can this be redirected to a string variable for comparison (and to suppress CLI output during testing?) end @testset "slicing" begin - xtal = deepcopy(sbmof1) + xtal = deepcopy(sbmof1) infer_bonds!(xtal, false) set_prop!(xtal.bonds, 1, 5, :bogus, "hi") ids = [1, 2, 5] @@ -147,298 +176,350 @@ end end @testset "rename" begin - xtal1 = deepcopy(sbmof1) - infer_bonds!(xtal1, true) - xtal2 = rename_xtal(xtal1, "renamed xtal") - - # test that crystal is unchanged - @test isapprox(xtal1.box, xtal2.box) - @test SimpleGraph(xtal1.bonds) == SimpleGraph(xtal2.bonds) - @test isapprox(xtal1.atoms, xtal2.atoms) - # test that name is changed - @test xtal2.name == "renamed xtal" + xtal1 = deepcopy(sbmof1) + infer_bonds!(xtal1, true) + xtal2 = rename_xtal(xtal1, "renamed xtal") + + # test that crystal is unchanged + @test isapprox(xtal1.box, xtal2.box) + @test SimpleGraph(xtal1.bonds) == SimpleGraph(xtal2.bonds) + @test isapprox(xtal1.atoms, xtal2.atoms) + # test that name is changed + @test xtal2.name == "renamed xtal" end @testset "cif reader" begin - # cif reader - xtal = deepcopy(test_structure2) - strip_numbers_from_atom_labels!(xtal) - @test xtal.name == "test_structure2.cif" - @test isapprox(xtal.box, Box(10.0, 20.0, 30.0, 90*π/180, 45*π/180, 120*π/180)) - @test xtal.atoms.n == 2 - @test isapprox(xtal.atoms, Atoms([:Ca, :O], Frac([0.2 0.6; 0.5 0.3; 0.7 0.1]))) - @test isapprox(xtal.charges, Charges([1.0, -1.0], Frac([0.2 0.6; 0.5 0.3; 0.7 0.1]))) - @test xtal.symmetry.is_p1 - @test Xtals.strip_number_from_label(:C) == :C - @test Xtals.strip_number_from_label(:C232) == :C - @test Xtals.strip_number_from_label(:Ca232) == :Ca - @test Xtals.strip_number_from_label(:Cad) == :Cad - - # test that incorrect file formats throw proper errors - @test_throws ErrorException Crystal("non_P1_no_symmetry.cif") - # test that a file with no atoms throws error - @test_throws ErrorException Crystal("no_atoms.cif") - - abs_path_xtal = Crystal(joinpath(rc[:paths][:crystals], "SBMOF-1.cif"), absolute_path=true) - @test sbmof1 ≈ abs_path_xtal + # cif reader + xtal = deepcopy(test_structure2) + strip_numbers_from_atom_labels!(xtal) + @test xtal.name == "test_structure2.cif" + @test isapprox( + xtal.box, + Box(10.0, 20.0, 30.0, 90 * π / 180, 45 * π / 180, 120 * π / 180) + ) + @test xtal.atoms.n == 2 + @test isapprox(xtal.atoms, Atoms([:Ca, :O], Frac([0.2 0.6; 0.5 0.3; 0.7 0.1]))) + @test isapprox(xtal.charges, Charges([1.0, -1.0], Frac([0.2 0.6; 0.5 0.3; 0.7 0.1]))) + @test xtal.symmetry.is_p1 + @test Xtals.strip_number_from_label(:C) == :C + @test Xtals.strip_number_from_label(:C232) == :C + @test Xtals.strip_number_from_label(:Ca232) == :Ca + @test Xtals.strip_number_from_label(:Cad) == :Cad + + # test that incorrect file formats throw proper errors + @test_throws ErrorException Crystal("non_P1_no_symmetry.cif") + # test that a file with no atoms throws error + @test_throws ErrorException Crystal("no_atoms.cif") + + abs_path_xtal = + Crystal(joinpath(rc[:paths][:crystals], "SBMOF-1.cif"); absolute_path=true) + @test sbmof1 ≈ abs_path_xtal end @testset "AtomsBase" begin - xtal = Crystal("SBMOF-1.cif") - pos = position(xtal) - @test length(pos) == 120 - @test isapprox(pos[1][2].val, 1.43954862, atol=1e-9) - @test ismissing(velocity(xtal)) - @test bounding_box(xtal)[1][1].val == xtal.box.a - @test periodicity(xtal) == [1,1,1] - @test n_dimensions(xtal) == 3 + xtal = Crystal("SBMOF-1.cif") + pos = position(xtal) + @test length(pos) == 120 + @test isapprox(pos[1][2].val, 1.43954862, atol=1e-9) + @test ismissing(velocity(xtal)) + @test bounding_box(xtal)[1][1].val == xtal.box.a + @test periodicity(xtal) == [1, 1, 1] + @test n_dimensions(xtal) == 3 end @testset "assign_charges" begin - xtal = Crystal("test_structure3.cif", include_zero_charges=false) - @test xtal.charges.n == 0 - strip_numbers_from_atom_labels!(xtal) - xtal2 = assign_charges(xtal, Dict(:Ca => -2.0, :O => 2.0)) - @test isapprox(xtal2.charges, Charges([-2.0, 2.0], Frac([0.2 0.6; 0.5 0.3; 0.7 0.1]))) - @test ! neutral(assign_charges(xtal, Dict(:Ca => 2.0, :O => 2.0), 100.0)) # not charge neutral - @test_throws ErrorException assign_charges(xtal, Dict(:Ca => 2.0, :O => 2.0)) # not charge neutral - @test empirical_formula(xtal) == Dict(:Ca => 1, :O => 1) - @test molecular_weight(xtal) ≈ 15.9994 + 40.078 - xtal3 = assign_charges(xtal2, Dict(:Ca => -2.0, :O => 2.0)) - @test isapprox(xtal2.charges, xtal3.charges) + xtal = Crystal("test_structure3.cif"; include_zero_charges=false) + @test xtal.charges.n == 0 + strip_numbers_from_atom_labels!(xtal) + xtal2 = assign_charges(xtal, Dict(:Ca => -2.0, :O => 2.0)) + @test isapprox(xtal2.charges, Charges([-2.0, 2.0], Frac([0.2 0.6; 0.5 0.3; 0.7 0.1]))) + @test !neutral(assign_charges(xtal, Dict(:Ca => 2.0, :O => 2.0), 100.0)) # not charge neutral + @test_throws ErrorException assign_charges(xtal, Dict(:Ca => 2.0, :O => 2.0)) # not charge neutral + @test empirical_formula(xtal) == Dict(:Ca => 1, :O => 1) + @test molecular_weight(xtal) ≈ 15.9994 + 40.078 + xtal3 = assign_charges(xtal2, Dict(:Ca => -2.0, :O => 2.0)) + @test isapprox(xtal2.charges, xtal3.charges) end @testset "overlap" begin - @test_throws ErrorException Crystal("test_structure2B.cif", check_overlap=true, check_neutrality=false) - @test_throws ErrorException Crystal("test_structure2B.cif", check_overlap=false, check_neutrality=true) - xtal = Crystal("test_structure2B.cif", check_overlap=false, check_neutrality=false, wrap_coords=true) - f = Frac([0.2 0.5 0.7; - 0.2 0.5 0.7; - 0.6 0.3 0.1; - 0.9 0.24 0.15]') - - @test isapprox(xtal.atoms, Atoms([:Ca, :Ca, :O, :C], f)) - @test isapprox(net_charge(xtal), 2.0) - @test empirical_formula(xtal) == Dict(:Ca => 2, :O => 1, :C => 1) - infer_bonds!(xtal, true) - @test_throws ErrorException Xtals.remove_duplicate_atoms_and_charges(xtal) - - @test_throws ErrorException Crystal("SBMOF-1_overlap.cif") - xtal1 = deepcopy(sbmof1_no_overlap) - @test isapprox(xtal1, sbmof1) - @test_throws ErrorException Crystal("SBMOF-1_overlap_diff_atom.cif", remove_duplicates=true) # duplicate but diff atom, so not repairable - - @test_throws ErrorException Crystal("example.mol") - @test_throws ErrorException Crystal("SBMOF-1.cif", infer_bonds=true) - infer_bonds!(xtal1, true) - @test_throws ErrorException replicate(xtal1, (1,1,1)) + @test_throws ErrorException Crystal( + "test_structure2B.cif", + check_overlap=true, + check_neutrality=false + ) + @test_throws ErrorException Crystal( + "test_structure2B.cif", + check_overlap=false, + check_neutrality=true + ) + xtal = Crystal( + "test_structure2B.cif"; + check_overlap=false, + check_neutrality=false, + wrap_coords=true + ) + f = Frac([ + 0.2 0.5 0.7 + 0.2 0.5 0.7 + 0.6 0.3 0.1 + 0.9 0.24 0.15 + ]') + + @test isapprox(xtal.atoms, Atoms([:Ca, :Ca, :O, :C], f)) + @test isapprox(net_charge(xtal), 2.0) + @test empirical_formula(xtal) == Dict(:Ca => 2, :O => 1, :C => 1) + infer_bonds!(xtal, true) + @test_throws ErrorException Xtals.remove_duplicate_atoms_and_charges(xtal) + + @test_throws ErrorException Crystal("SBMOF-1_overlap.cif") + xtal1 = deepcopy(sbmof1_no_overlap) + @test isapprox(xtal1, sbmof1) + @test_throws ErrorException Crystal( + "SBMOF-1_overlap_diff_atom.cif", + remove_duplicates=true + ) # duplicate but diff atom, so not repairable + + @test_throws ErrorException Crystal("example.mol") + @test_throws ErrorException Crystal("SBMOF-1.cif", infer_bonds=true) + infer_bonds!(xtal1, true) + @test_throws ErrorException replicate(xtal1, (1, 1, 1)) end @testset "xyz writer" begin - xtal = deepcopy(sbmof1) - xtal_xyz_filename = replace(replace(xtal.name, ".cif" => ""), ".cssr" => "") * ".xyz" - write_xyz(xtal) - atoms_read = read_xyz(xtal_xyz_filename) # Atoms{Cart} - atoms_read_f = Frac(atoms_read, xtal.box) # Atoms{Frac} - @test isapprox(atoms_read_f, xtal.atoms, atol=0.001) - write_xyz(xtal, center_at_origin=true, comment="blah") #center coords - atoms_read = read_xyz(xtal_xyz_filename) # Atoms{Cart} - atoms_read_f = Frac(atoms_read, xtal.box) # Atoms{Frac} - @test isapprox(atoms_read_f.coords.xf, xtal.atoms.coords.xf .- [0.5, 0.5, 0.5], atol=0.001) - rm(xtal_xyz_filename) # clean up + xtal = deepcopy(sbmof1) + xtal_xyz_filename = replace(replace(xtal.name, ".cif" => ""), ".cssr" => "") * ".xyz" + write_xyz(xtal) + atoms_read = read_xyz(xtal_xyz_filename) # Atoms{Cart} + atoms_read_f = Frac(atoms_read, xtal.box) # Atoms{Frac} + @test isapprox(atoms_read_f, xtal.atoms, atol=0.001) + write_xyz(xtal; center_at_origin=true, comment="blah") #center coords + atoms_read = read_xyz(xtal_xyz_filename) # Atoms{Cart} + atoms_read_f = Frac(atoms_read, xtal.box) # Atoms{Frac} + @test isapprox( + atoms_read_f.coords.xf, + xtal.atoms.coords.xf .- [0.5, 0.5, 0.5], + atol=0.001 + ) + rm(xtal_xyz_filename) # clean up end @testset "cif writer" begin - crystal = deepcopy(sbmof1) - rewrite_filename = "rewritten.cif" - write_cif(crystal, joinpath(rc[:paths][:crystals], rewrite_filename)) - crystal_reloaded = Crystal(rewrite_filename) - strip_numbers_from_atom_labels!(crystal_reloaded) # TODO Arthur remove this necessity from write_cif - @test isapprox(crystal, crystal_reloaded, atol=0.0001) - write_cif(crystal, joinpath(rc[:paths][:crystals], rewrite_filename), fractional_coords=false) # cartesian - crystal_reloaded = Crystal(rewrite_filename) - strip_numbers_from_atom_labels!(crystal_reloaded) # TODO Arthur remove this necessity from write_cif - @test isapprox(crystal, crystal_reloaded, atol=0.0001) - rewrite_filename = "rewritten.cif" - write_cif(crystal, joinpath(rc[:paths][:crystals], rewrite_filename), number_atoms=false) - crystal_reloaded = Crystal(rewrite_filename) - @test isapprox(crystal, crystal_reloaded, atol=0.0001) - crystal = Crystal("ATIBOU01_clean.cssr", include_zero_charges=false) - @test_throws ErrorException write_cif(crystal) - crystal = assign_charges(crystal, Dict(:Mn => 2.0, :O => -2.0, :H => 0.0, :C => 0.0), Inf) - write_cif(crystal) - @test isfile("ATIBOU01_clean.cssr.cif") + crystal = deepcopy(sbmof1) + rewrite_filename = "rewritten.cif" + write_cif(crystal, joinpath(rc[:paths][:crystals], rewrite_filename)) + crystal_reloaded = Crystal(rewrite_filename) + strip_numbers_from_atom_labels!(crystal_reloaded) # TODO Arthur remove this necessity from write_cif + @test isapprox(crystal, crystal_reloaded, atol=0.0001) + write_cif( + crystal, + joinpath(rc[:paths][:crystals], rewrite_filename); + fractional_coords=false + ) # cartesian + crystal_reloaded = Crystal(rewrite_filename) + strip_numbers_from_atom_labels!(crystal_reloaded) # TODO Arthur remove this necessity from write_cif + @test isapprox(crystal, crystal_reloaded, atol=0.0001) + rewrite_filename = "rewritten.cif" + write_cif( + crystal, + joinpath(rc[:paths][:crystals], rewrite_filename); + number_atoms=false + ) + crystal_reloaded = Crystal(rewrite_filename) + @test isapprox(crystal, crystal_reloaded, atol=0.0001) + crystal = Crystal("ATIBOU01_clean.cssr"; include_zero_charges=false) + @test_throws ErrorException write_cif(crystal) + crystal = + assign_charges(crystal, Dict(:Mn => 2.0, :O => -2.0, :H => 0.0, :C => 0.0), Inf) + write_cif(crystal) + @test isfile("ATIBOU01_clean.cssr.cif") end @testset "strip numbers from labels" begin - hk = Crystal("HKUST-1_low_symm.cif", remove_duplicates=true) - strip_numbers_from_atom_labels!(hk) - hk_p1 = Crystal("HKUST-1_P1.cif") # from avogadro - strip_numbers_from_atom_labels!(hk_p1) - @test equal_multisets(hk.atoms, hk_p1.atoms, hk.box) + hk = Crystal("HKUST-1_low_symm.cif"; remove_duplicates=true) + strip_numbers_from_atom_labels!(hk) + hk_p1 = Crystal("HKUST-1_P1.cif") # from avogadro + strip_numbers_from_atom_labels!(hk_p1) + @test equal_multisets(hk.atoms, hk_p1.atoms, hk.box) end @testset "cssr read/write" begin - xtal_cssr = Crystal("test_structure2.cssr") - strip_numbers_from_atom_labels!(xtal_cssr) - xtal_cif = deepcopy(test_structure2) - strip_numbers_from_atom_labels!(xtal_cif) - @test isapprox(xtal_cif, xtal_cssr) - - da_xtalname = "ATIBOU01_clean" - xtal_cif = Crystal(da_xtalname * ".cif", include_zero_charges=true) - write_cssr(xtal_cif, quiet=true) - write_cssr(xtal_cif, joinpath(rc[:paths][:crystals], da_xtalname * ".cssr")) - xtal_cssr = Crystal(da_xtalname * ".cssr", include_zero_charges=true) - @test isapprox(xtal_cssr, xtal_cif, atol=0.0001) + xtal_cssr = Crystal("test_structure2.cssr") + strip_numbers_from_atom_labels!(xtal_cssr) + xtal_cif = deepcopy(test_structure2) + strip_numbers_from_atom_labels!(xtal_cif) + @test isapprox(xtal_cif, xtal_cssr) + + da_xtalname = "ATIBOU01_clean" + xtal_cif = Crystal(da_xtalname * ".cif"; include_zero_charges=true) + write_cssr(xtal_cif; quiet=true) + write_cssr(xtal_cif, joinpath(rc[:paths][:crystals], da_xtalname * ".cssr")) + xtal_cssr = Crystal(da_xtalname * ".cssr"; include_zero_charges=true) + @test isapprox(xtal_cssr, xtal_cif, atol=0.0001) end @testset "replication" begin - sbmof = deepcopy(sbmof1) - strip_numbers_from_atom_labels!(sbmof) - replicated_sbmof = replicate(sbmof, (1, 1, 1)) - @test isapprox(sbmof, replicated_sbmof) - @test isapprox(2 * 2 * molecular_weight(sbmof), molecular_weight(replicate(sbmof, (2, 2, 1)))) - @test isapprox(2 * sbmof.box.a, replicate(sbmof, (2, 3, 4)).box.a) - - repfactors = Xtals.replication_factors(sbmof.box, 14.0) - replicated_sbmof = replicate(sbmof, repfactors) - @test Xtals.replication_factors(replicated_sbmof.box, 14.0) == (1, 1, 1) - @test isapprox(sbmof.atoms.coords.xf[:, 1] ./ repfactors, replicated_sbmof.atoms.coords.xf[:, 1]) - @test isapprox(replicated_sbmof.box.reciprocal_lattice, 2 * π * inv(replicated_sbmof.box.f_to_c)) - @test empirical_formula(sbmof) == empirical_formula(replicated_sbmof) - @test isapprox(crystal_density(sbmof), crystal_density(replicated_sbmof), atol=1e-7) - - xtal = deepcopy(sbmof1) - rbox = replicate(xtal.box, (2, 3, 4)) - @test rbox.Ω ≈ xtal.box.Ω * 2 * 3 * 4 - @test all(rbox.c_to_f * xtal.box.f_to_c * [1.0, 1.0, 1.0] .≈ [1/2, 1/3, 1/4]) + sbmof = deepcopy(sbmof1) + strip_numbers_from_atom_labels!(sbmof) + replicated_sbmof = replicate(sbmof, (1, 1, 1)) + @test isapprox(sbmof, replicated_sbmof) + @test isapprox( + 2 * 2 * molecular_weight(sbmof), + molecular_weight(replicate(sbmof, (2, 2, 1))) + ) + @test isapprox(2 * sbmof.box.a, replicate(sbmof, (2, 3, 4)).box.a) + + repfactors = Xtals.replication_factors(sbmof.box, 14.0) + replicated_sbmof = replicate(sbmof, repfactors) + @test Xtals.replication_factors(replicated_sbmof.box, 14.0) == (1, 1, 1) + @test isapprox( + sbmof.atoms.coords.xf[:, 1] ./ repfactors, + replicated_sbmof.atoms.coords.xf[:, 1] + ) + @test isapprox( + replicated_sbmof.box.reciprocal_lattice, + 2 * π * inv(replicated_sbmof.box.f_to_c) + ) + @test empirical_formula(sbmof) == empirical_formula(replicated_sbmof) + @test isapprox(crystal_density(sbmof), crystal_density(replicated_sbmof), atol=1e-7) + + xtal = deepcopy(sbmof1) + rbox = replicate(xtal.box, (2, 3, 4)) + @test rbox.Ω ≈ xtal.box.Ω * 2 * 3 * 4 + @test all(rbox.c_to_f * xtal.box.f_to_c * [1.0, 1.0, 1.0] .≈ [1 / 2, 1 / 3, 1 / 4]) end @testset "vtk writer" begin - xtal = deepcopy(sbmof1) - @test Xtals.vtk_filename(xtal) == "SBMOF-1.vtk" - xtal = Crystal("ATIBOU01_clean.cssr") - @test Xtals.vtk_filename(xtal) == "ATIBOU01_clean.vtk" + xtal = deepcopy(sbmof1) + @test Xtals.vtk_filename(xtal) == "SBMOF-1.vtk" + xtal = Crystal("ATIBOU01_clean.cssr") + @test Xtals.vtk_filename(xtal) == "ATIBOU01_clean.vtk" end @testset "addition" begin - # test crystal addition - c1 = Crystal("crystal 1", unit_cube(), Atoms([:a, :b], - Frac([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) - ), - Charges([0.1, 0.2], - Frac([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) - ) - ) - c2 = Crystal("crystal 2", unit_cube(), Atoms([:c, :d], - Frac([7.0 10.0; - 8.0 11.0; - 9.0 12.0] - ), - ), - Charges([0.3, 0.4], - Frac([7.0 10.0; - 8.0 11.0; - 9.0 12.0] - ) - ) - ) - c = c1 + c2 - @test_throws AssertionError c1 + sbmof1 # only allow crystals with same box - @test isapprox(c1.box, c.box) - @test isapprox(c2.box, c.box) - @test isapprox(c1.atoms + c2.atoms, c.atoms) - @test isapprox(c1.charges + c2.charges, c.charges) - @test isapprox(c[1:2], c1) # indexing test - @test isapprox(c[3:4], c2) # indexing test - - # test overlap crystal addition - c_overlap = +(c1, c2, c; check_overlap=false) - @test isapprox(c_overlap.box, c.box) - @test isapprox(c_overlap.atoms, c1.atoms + c2.atoms + c.atoms) - @test isapprox(c_overlap.charges, c1.charges + c2.charges + c.charges) - @test_throws AssertionError +(c1, c2, c) # by default, do not let crystal additions result in overlap - - # bond addition - xtal = deepcopy(sbmof1) - xtal2 = deepcopy(xtal) - infer_bonds!(xtal, true) - infer_bonds!(xtal2, true) - xtal3 = +(xtal, xtal2; check_overlap=false) - @test ne(xtal3.bonds) == ne(xtal.bonds) + ne(xtal2.bonds) + # test crystal addition + c1 = Crystal( + "crystal 1", + unit_cube(), + Atoms([:a, :b], Frac([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ])), + Charges([0.1, 0.2], Frac([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ])) + ) + c2 = Crystal( + "crystal 2", + unit_cube(), + Atoms([:c, :d], Frac([ + 7.0 10.0 + 8.0 11.0 + 9.0 12.0 + ])), + Charges([0.3, 0.4], Frac([ + 7.0 10.0 + 8.0 11.0 + 9.0 12.0 + ])) + ) + c = c1 + c2 + @test_throws AssertionError c1 + sbmof1 # only allow crystals with same box + @test isapprox(c1.box, c.box) + @test isapprox(c2.box, c.box) + @test isapprox(c1.atoms + c2.atoms, c.atoms) + @test isapprox(c1.charges + c2.charges, c.charges) + @test isapprox(c[1:2], c1) # indexing test + @test isapprox(c[3:4], c2) # indexing test + + # test overlap crystal addition + c_overlap = +(c1, c2, c; check_overlap=false) + @test isapprox(c_overlap.box, c.box) + @test isapprox(c_overlap.atoms, c1.atoms + c2.atoms + c.atoms) + @test isapprox(c_overlap.charges, c1.charges + c2.charges + c.charges) + @test_throws AssertionError +(c1, c2, c) # by default, do not let crystal additions result in overlap + + # bond addition + xtal = deepcopy(sbmof1) + xtal2 = deepcopy(xtal) + infer_bonds!(xtal, true) + infer_bonds!(xtal2, true) + xtal3 = +(xtal, xtal2; check_overlap=false) + @test ne(xtal3.bonds) == ne(xtal.bonds) + ne(xtal2.bonds) end @testset "indexing" begin - sbmof1_sub = sbmof1[10:15] - @test sbmof1_sub.atoms.n == 6 - @test sbmof1_sub.charges.n == 6 - @test sbmof1_sub.atoms.species == sbmof1.atoms.species[10:15] - @test sbmof1_sub.atoms.coords.xf == sbmof1.atoms.coords.xf[:, 10:15] - unequal_n_q = Crystal("symmetry_test_structure_charges.cif", include_zero_charges=false) - @test_throws ErrorException lastindex(unequal_n_q) - @test_throws ErrorException unequal_n_q[[1]] - - @test !isapprox(unequal_n_q, symmetry_test_structure) + sbmof1_sub = sbmof1[10:15] + @test sbmof1_sub.atoms.n == 6 + @test sbmof1_sub.charges.n == 6 + @test sbmof1_sub.atoms.species == sbmof1.atoms.species[10:15] + @test sbmof1_sub.atoms.coords.xf == sbmof1.atoms.coords.xf[:, 10:15] + unequal_n_q = Crystal("symmetry_test_structure_charges.cif"; include_zero_charges=false) + @test_throws ErrorException lastindex(unequal_n_q) + @test_throws ErrorException unequal_n_q[[1]] + + @test !isapprox(unequal_n_q, symmetry_test_structure) end @testset "charges" begin - # including zero charges or not, when reading in a .cif `include_zero_charges` flag to Crystal constructor - frame1 = Crystal("ATIBOU01_clean.cif", include_zero_charges=false) # has four zero charges in it - @test ! any(frame1.charges.q .== 0.0) - @test frame1.charges.n == frame1.atoms.n - 4 # 4 charges are zero - frame2 = Crystal("ATIBOU01_clean.cif") - @test frame2.charges.n == frame2.atoms.n - @test isapprox(frame2.charges.coords, frame2.atoms.coords) - - crystal = replicate(Crystal("CAXVII_clean.cif"), (2, 2, 2)) - @test isapprox(sum(crystal.charges.q), 0.0, atol=0.001) - - # high-precision charges in write_cif - crystal = Crystal("ADARAA_clean.cif") - @test all(crystal.charges.q .!= 0.0) # has charges... - @test neutral(crystal, Xtals.NET_CHARGE_TOL) # is neutral... - write_cif(crystal, joinpath(rc[:paths][:crystals], "high_precision_charges_test.cif")) # should write charges with more decimals - crystal_reloaded = Crystal("high_precision_charges_test.cif") - @test isapprox(crystal.charges.q, crystal_reloaded.charges.q, atol=1e-8) # would not hv this precision if didn't write high-precision charges for this one. - crystal_not_neutral = crystal.charges.q .= round.(crystal.charges.q, digits=6) # this is what charges would be if stored 6 decimals in .cif and threw away the rest - write_cif(crystal, joinpath(rc[:paths][:crystals], "high_precision_charges_test_not_neutral.cif")) - @test ! neutral(Crystal("high_precision_charges_test_not_neutral.cif", check_neutrality=false), Xtals.NET_CHARGE_TOL) # not gonna be neutral - @test_throws ErrorException Crystal("high_precision_charges_test_not_neutral.cif") # should throw error b/c not charge neutral + # including zero charges or not, when reading in a .cif `include_zero_charges` flag to Crystal constructor + frame1 = Crystal("ATIBOU01_clean.cif"; include_zero_charges=false) # has four zero charges in it + @test !any(frame1.charges.q .== 0.0) + @test frame1.charges.n == frame1.atoms.n - 4 # 4 charges are zero + frame2 = Crystal("ATIBOU01_clean.cif") + @test frame2.charges.n == frame2.atoms.n + @test isapprox(frame2.charges.coords, frame2.atoms.coords) + + crystal = replicate(Crystal("CAXVII_clean.cif"), (2, 2, 2)) + @test isapprox(sum(crystal.charges.q), 0.0, atol=0.001) + + # high-precision charges in write_cif + crystal = Crystal("ADARAA_clean.cif") + @test all(crystal.charges.q .!= 0.0) # has charges... + @test neutral(crystal, Xtals.NET_CHARGE_TOL) # is neutral... + write_cif(crystal, joinpath(rc[:paths][:crystals], "high_precision_charges_test.cif")) # should write charges with more decimals + crystal_reloaded = Crystal("high_precision_charges_test.cif") + @test isapprox(crystal.charges.q, crystal_reloaded.charges.q, atol=1e-8) # would not hv this precision if didn't write high-precision charges for this one. + crystal_not_neutral = crystal.charges.q .= round.(crystal.charges.q, digits=6) # this is what charges would be if stored 6 decimals in .cif and threw away the rest + write_cif( + crystal, + joinpath(rc[:paths][:crystals], "high_precision_charges_test_not_neutral.cif") + ) + @test !neutral( + Crystal("high_precision_charges_test_not_neutral.cif"; check_neutrality=false), + Xtals.NET_CHARGE_TOL + ) # not gonna be neutral + @test_throws ErrorException Crystal("high_precision_charges_test_not_neutral.cif") # should throw error b/c not charge neutral end @testset "species_from_col" begin - # tests for species_from_col: - # 1. Make sure a .cif loads a Crystal with bad labels as default; - # fixed labels with species_column="_atom_site_type_symbol". - # 2. Make sure a .cif w/ good labels loads to an identical Crystal with - # and without species_column="_atom_site_type_symbol". - # 3. Make sure a bogus species_column throws an exception. - xtal = Crystal("SBMOF-1_bad_labels.cif", species_col=["_atom_site_label"]) - @test xtal.atoms.species[1] == :C1A - - xtal = Crystal("SBMOF-1_bad_labels.cif", species_col=["_atom_site_type_symbol"]) - @test xtal.atoms.species[1] == :C - - xtal1 = deepcopy(sbmof1) - xtal2 = deepcopy(sbmof1) - @test xtal1.atoms.species == xtal2.atoms.species - - @test_throws KeyError Crystal("SBMOF-1.cif", species_col=["bogus_column"]) + # tests for species_from_col: + # 1. Make sure a .cif loads a Crystal with bad labels as default; + # fixed labels with species_column="_atom_site_type_symbol". + # 2. Make sure a .cif w/ good labels loads to an identical Crystal with + # and without species_column="_atom_site_type_symbol". + # 3. Make sure a bogus species_column throws an exception. + xtal = Crystal("SBMOF-1_bad_labels.cif"; species_col=["_atom_site_label"]) + @test xtal.atoms.species[1] == :C1A + + xtal = Crystal("SBMOF-1_bad_labels.cif"; species_col=["_atom_site_type_symbol"]) + @test xtal.atoms.species[1] == :C + + xtal1 = deepcopy(sbmof1) + xtal2 = deepcopy(sbmof1) + @test xtal1.atoms.species == xtal2.atoms.species + + @test_throws KeyError Crystal("SBMOF-1.cif", species_col=["bogus_column"]) end @testset "other" begin - xtal = Crystal("SBMOF-1.cif", include_zero_charges=false) - @test ! has_charges(xtal) - @test isapprox(xtal.box.reciprocal_lattice, 2 * π * inv(xtal.box.f_to_c)) - @test xtal.box.Ω ≈ det(xtal.box.f_to_c) # sneak in crystal test - @test isapprox(crystal_density(xtal), 1570.4, atol=0.5) # kg/m3 - end + xtal = Crystal("SBMOF-1.cif"; include_zero_charges=false) + @test !has_charges(xtal) + @test isapprox(xtal.box.reciprocal_lattice, 2 * π * inv(xtal.box.f_to_c)) + @test xtal.box.Ω ≈ det(xtal.box.f_to_c) # sneak in crystal test + @test isapprox(crystal_density(xtal), 1570.4, atol=0.5) # kg/m3 +end end diff --git a/test/distance.jl b/test/distance.jl index 1ac833f..426420d 100644 --- a/test/distance.jl +++ b/test/distance.jl @@ -14,32 +14,34 @@ using LinearAlgebra @test isapprox(dxf, [-0.3, -0.1, 0.1]) # distance (fractional) - f = Frac([0.2 0.4; - 0.1 0.8; - 0.8 0.6] - ) + f = Frac([ + 0.2 0.4 + 0.1 0.8 + 0.8 0.6 + ]) box = unit_cube() - @test isapprox(distance(f, box, 1, 2, false), norm(f.xf[:, 1] - f.xf[:, 2])) + @test isapprox(distance(f, box, 1, 2, false), norm(f.xf[:, 1] - f.xf[:, 2])) @test isapprox(distance(f, box, 2, 2, false), 0.0) @test isapprox(distance(f, box, 2, 1, false), distance(f, box, 1, 2, false)) - @test isapprox(distance(f, box, 1, 2, true), norm(f.xf[:, 1] - [0.4, -0.2, 0.6])) + @test isapprox(distance(f, box, 1, 2, true), norm(f.xf[:, 1] - [0.4, -0.2, 0.6])) box = Box(1.0, 10.0, 100.0) - @test isapprox(distance(f, box, 1, 2, false), norm([0.2, 1.0, 80.0] - [0.4, 8.0, 60.0])) - @test isapprox(distance(f, box, 2, 1, true), norm([0.2, 1.0, 80.0] - [0.4, -2.0, 60.0])) + @test isapprox(distance(f, box, 1, 2, false), norm([0.2, 1.0, 80.0] - [0.4, 8.0, 60.0])) + @test isapprox(distance(f, box, 2, 1, true), norm([0.2, 1.0, 80.0] - [0.4, -2.0, 60.0])) # distance (Cartesian) - c = Cart([0.2 0.4; - 0.1 0.8; - 0.8 0.6] - ) + c = Cart([ + 0.2 0.4 + 0.1 0.8 + 0.8 0.6 + ]) box = unit_cube() - @test isapprox(distance(c, box, 1, 2, false), norm(c.x[:, 1] - c.x[:, 2])) + @test isapprox(distance(c, box, 1, 2, false), norm(c.x[:, 1] - c.x[:, 2])) @test isapprox(distance(c, box, 2, 2, false), 0.0) @test isapprox(distance(c, box, 2, 1, false), distance(c, box, 1, 2, false)) - @test isapprox(distance(c, box, 1, 2, true), norm(c.x[:, 1] - [0.4, -0.2, 0.6])) + @test isapprox(distance(c, box, 1, 2, true), norm(c.x[:, 1] - [0.4, -0.2, 0.6])) box = Box(10.0, 1.0, 100.0) @test isapprox(distance(c, box, 1, 2, false), norm(c.x[:, 1] - c.x[:, 2])) - @test isapprox(distance(c, box, 2, 1, true), norm([0.2, 0.1, 0.8] - [0.4, -0.2, 0.6])) + @test isapprox(distance(c, box, 2, 1, true), norm([0.2, 0.1, 0.8] - [0.4, -0.2, 0.6])) # distance (Tests from Avogadro measurement tool) crystal = Crystal("distance_test.cif") @@ -50,32 +52,34 @@ using LinearAlgebra # overlap box = unit_cube() - f = Frac([0.2 0.4 0.6 0.401; - 0.1 0.8 0.7 0.799; - 0.8 0.6 0.5 0.602] - ) - o_flag, o_ids = overlap(f, box, true, tol=0.01) + f = Frac([ + 0.2 0.4 0.6 0.401 + 0.1 0.8 0.7 0.799 + 0.8 0.6 0.5 0.602 + ]) + o_flag, o_ids = overlap(f, box, true; tol=0.01) @test o_flag @test o_ids == [(2, 4)] - f = Frac([0.2 0.4 0.6 0.2; - 0.1 0.8 0.7 0.1; - 0.99 0.6 0.5 0.01] - ) - o_flag, o_ids = overlap(f, box, true, tol=0.03) + f = Frac([ + 0.2 0.4 0.6 0.2 + 0.1 0.8 0.7 0.1 + 0.99 0.6 0.5 0.01 + ]) + o_flag, o_ids = overlap(f, box, true; tol=0.03) @test o_flag @test o_ids == [(1, 4)] - o_flag, o_ids = overlap(f, box, false, tol=0.03) + o_flag, o_ids = overlap(f, box, false; tol=0.03) @test !o_flag @test o_ids == [] xtal = Crystal("IRMOF-1.cif") @test !overlap(xtal, true)[1] - xtal = Crystal("IRMOF-1_overlap.cif", check_overlap=false) + xtal = Crystal("IRMOF-1_overlap.cif"; check_overlap=false) @test overlap(xtal, true)[1] # test distance function (via Avogadro) - crystal = Crystal("simple_test.cif", check_overlap=false) + crystal = Crystal("simple_test.cif"; check_overlap=false) @test distance(crystal.atoms, crystal.box, 1, 1, true) == 0.0 @test isapprox(distance(crystal.atoms, crystal.box, 2, 5, true), 4.059, atol=0.001) @test isapprox(distance(crystal.atoms, crystal.box, 2, 5, false), 4.059, atol=0.001) @@ -122,22 +126,22 @@ using LinearAlgebra # vector-column distances ### - function check_vector_indices(coords, box, Is, Js, apply_pbc = true) - d = distance(coords, box, Is, Js, apply_pbc) + function check_vector_indices(coords, box, Is, Js, apply_pbc=true) + d = distance(coords, box, Is, Js, apply_pbc) - # force using distance(::Coords, ::Box, ::Integer, ::Integer) - d_elementwise = distance.(Ref(coords), Ref(box), Is, Js, apply_pbc) + # force using distance(::Coords, ::Box, ::Integer, ::Integer) + d_elementwise = distance.(Ref(coords), Ref(box), Is, Js, apply_pbc) - # CartesianIndex is implicity a 1-vector - if Is isa CartesianIndices - d_elementwise = first.(d_elementwise) - end + # CartesianIndex is implicity a 1-vector + if Is isa CartesianIndices + d_elementwise = first.(d_elementwise) + end - @test d ≈ d_elementwise + @test d ≈ d_elementwise end # different types of indices - Is = (1:5, rand(1:5, 5), CartesianIndices((1:5))) + Is = (1:5, rand(1:5, 5), CartesianIndices((1:5))) Js = (1:5, rand(1:5, 5), CartesianIndices((1:5))) for (I, J) in zip(Is, Js), pbc in (true, false) diff --git a/test/matter.jl b/test/matter.jl index f5f35fe..a4a8fd2 100644 --- a/test/matter.jl +++ b/test/matter.jl @@ -4,25 +4,29 @@ using Xtals using Test @testset "Matter Tests" begin - f1 = Frac([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) + f1 = Frac([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ]) @test length(f1) == 2 - f2 = Frac([7.0 10.0; - 8.0 11.0; - 9.0 12.0] - ) - - c1 = Cart([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) + f2 = Frac([ + 7.0 10.0 + 8.0 11.0 + 9.0 12.0 + ]) + + c1 = Cart([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ]) @test length(c1) == 2 - c2 = Cart([7.0 10.0; - 8.0 11.0; - 9.0 12.0] - ) + c2 = Cart([ + 7.0 10.0 + 8.0 11.0 + 9.0 12.0 + ]) @test isapprox(c2[1], Cart([7.0; 8.0; 9.0])) # test getindex(::Cart,::Int) @test lastindex(c2) == 2 c2[2] = [0.0, 1.0, 2.0] @@ -38,12 +42,12 @@ using Test a1_c = Atoms(s1, c1) a2_c = Atoms(s2, c2) - @test ! isapprox(a1_c, a2_c) - @test ! isapprox(a1_f, a2_f) - @test ! isapprox(Atoms(s1, f1), Atoms(s1, f2)) - @test ! isapprox(Atoms(s1, f1), Atoms(s2, f1)) - @test ! isapprox(Atoms(s1, c1), Atoms(s1, c2)) - @test ! isapprox(Atoms(s1, c1), Atoms(s2, c1)) + @test !isapprox(a1_c, a2_c) + @test !isapprox(a1_f, a2_f) + @test !isapprox(Atoms(s1, f1), Atoms(s1, f2)) + @test !isapprox(Atoms(s1, f1), Atoms(s2, f1)) + @test !isapprox(Atoms(s1, c1), Atoms(s1, c2)) + @test !isapprox(Atoms(s1, c1), Atoms(s2, c1)) a3_c = a1_c + a2_c @test a3_c.n == a1_c.n + a2_c.n @@ -65,17 +69,17 @@ using Test q1_c = Charges(q_vals_1, c1) q2_c = Charges(q_vals_2, c2) - @test ! neutral(q1_c) + @test !neutral(q1_c) q1_f = Charges(q_vals_1, f1) q2_f = Charges(q_vals_2, f2) - @test ! isapprox(q1_c, q2_c) - @test ! isapprox(q1_f, q2_f) - @test ! isapprox(Charges(q_vals_1, f1), Charges(q_vals_1, f2)) - @test ! isapprox(Charges(q_vals_1, f1), Charges(q_vals_2, f1)) - @test ! isapprox(Charges(q_vals_1, c1), Charges(q_vals_1, c2)) - @test ! isapprox(Charges(q_vals_1, c1), Charges(q_vals_2, c1)) + @test !isapprox(q1_c, q2_c) + @test !isapprox(q1_f, q2_f) + @test !isapprox(Charges(q_vals_1, f1), Charges(q_vals_1, f2)) + @test !isapprox(Charges(q_vals_1, f1), Charges(q_vals_2, f1)) + @test !isapprox(Charges(q_vals_1, c1), Charges(q_vals_1, c2)) + @test !isapprox(Charges(q_vals_1, c1), Charges(q_vals_2, c1)) q3_c = q1_c + q2_c @test q3_c.n == q1_c.n + q2_c.n @@ -92,14 +96,16 @@ using Test @test q3_f.coords.xf[:, 3:4] == f2.xf # wrap - f = Frac([0.1 -0.1; - 2.2 0.56; - -0.4 1.2] - ) - fw = Frac([0.1 0.9; - 0.2 0.56; - 0.6 0.2] - ) + f = Frac([ + 0.1 -0.1 + 2.2 0.56 + -0.4 1.2 + ]) + fw = Frac([ + 0.1 0.9 + 0.2 0.56 + 0.6 0.2 + ]) wrap!(f) @test isapprox(f, fw) @@ -107,7 +113,7 @@ using Test @test neutral(q) @test isapprox(net_charge(q), 0.0) q.q[2] = 0.0 - @test ! neutral(q) + @test !neutral(q) @test isapprox(net_charge(q), -0.2) q = Charges{Frac}(0) @test net_charge(q) == 0.0 @@ -143,29 +149,30 @@ using Test @test length(charges[6:10].q) == 5 == charges[6:10].n # translate_by - f = Frac([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) + f = Frac([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ]) dxf = Frac([1.0, 3.0, -1.0]) translate_by!(f, dxf) - @test isapprox(f, Frac([2.0 5.0; - 5.0 8.0; - 2.0 5.0] - ) - ) - - c = Cart([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) + @test isapprox(f, Frac([ + 2.0 5.0 + 5.0 8.0 + 2.0 5.0 + ])) + + c = Cart([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ]) dx = Cart([1.0, 3.0, -1.0]) translate_by!(c, dx) - @test isapprox(c, Cart([2.0 5.0; - 5.0 8.0; - 2.0 5.0] - ) - ) - + @test isapprox(c, Cart([ + 2.0 5.0 + 5.0 8.0 + 2.0 5.0 + ])) end end diff --git a/test/misc.jl b/test/misc.jl index 44a8579..2b211f7 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -15,10 +15,11 @@ import IOCapture.capture @test isapprox(am[:Co], 58.9332, atol=0.001) test_xyz_filename = "atoms_test" - c = Cart([1.0 4.0; - 2.0 5.0; - 3.0 6.0] - ) + c = Cart([ + 1.0 4.0 + 2.0 5.0 + 3.0 6.0 + ]) s = [:C, :H] atoms = Atoms(s, c) write_xyz(atoms, test_xyz_filename) @@ -26,18 +27,20 @@ import IOCapture.capture @test isapprox(atoms, atoms_read) rm(test_xyz_filename * ".xyz") - @test rc[:cpk_colors][:Li] == (204,128,255) + @test rc[:cpk_colors][:Li] == (204, 128, 255) atoms, bonds, bond_types = read_mol("data/example.mol") - @test (atoms.species[1] == :O) && (atoms.species[end] == :H) & (length(atoms.species) == 31) && (atoms.n == 31) + @test (atoms.species[1] == :O) && + (atoms.species[end] == :H) & (length(atoms.species) == 31) && + (atoms.n == 31) @test ne(bonds) == 32 @test nv(bonds) == 31 @test (bond_types[1] == 1) && (bond_types[3] == 2) @test neighbors(bonds, 1) == [2, 5] - xtal = Crystal("SBMOF-1.cif", infer_bonds=true, periodic_boundaries=true) + xtal = Crystal("SBMOF-1.cif"; infer_bonds=true, periodic_boundaries=true) mol2_temp = tempname() - write_mol2(xtal, filename=mol2_temp) + write_mol2(xtal; filename=mol2_temp) @test isfile(mol2_temp) write_mol2(xtal) @test isfile("SBMOF-1.mol2") diff --git a/test/paths.jl b/test/paths.jl index 96e4c63..9f6364b 100644 --- a/test/paths.jl +++ b/test/paths.jl @@ -6,7 +6,7 @@ using IOCapture, Test, Xtals @test rc[:atomic_masses][:He] == 4.0026 oldpath = rc[:paths][:data] newpath = joinpath(pwd(), "other_data") - set_paths(newpath, print_paths=true) + set_paths(newpath; print_paths=true) @test rc[:paths][:crystals] == joinpath(newpath, "crystals") rc[:paths][:crystals] = joinpath(newpath, "other_crystals") xtal = Crystal("SBMOF-1.cif") diff --git a/test/runtests.jl b/test/runtests.jl index d1be604..f690b11 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ testfiles = [ "box.jl" "assert_p1_symmetry.jl" "paths.jl" - ] +] @assert (VERSION.major == 1) && (VERSION.minor ≥ 6) "Minimum Julia version not met." @@ -15,10 +15,10 @@ using Documenter, IOCapture, Logging, Test, Xtals Xtals.banner() -for testfile ∈ testfiles +for testfile in testfiles @info "Running test/$testfile" with_logger(NullLogger()) do - include(testfile) + return include(testfile) end end