diff --git a/.gitignore b/.gitignore index 285c8fcb1..c314cf993 100644 --- a/.gitignore +++ b/.gitignore @@ -8,4 +8,6 @@ *.jld2 docs/build/ docs/site/ -*checkpoint.ipynb \ No newline at end of file +*checkpoint.ipynb +.vscode +Manifest.toml diff --git a/.travis.yml b/.travis.yml index 9ea90275c..40531c6a1 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,9 +1,9 @@ language: julia -coveralls: true +codecov: true julia: - - 1.5 + - 1.6 addons: apt: @@ -13,20 +13,13 @@ addons: sudo: enabled script: - - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); Pkg.test("PorousMaterials")' + - julia --project -e 'ENV["JULIA_PKG_SERVER"] = ""; import Pkg; Pkg.update()' + - julia --project --check-bounds=yes -e 'import Pkg; Pkg.build(); Pkg.test("PorousMaterials")' after_success: - - julia --project -e 'import Pkg; Pkg.add("Documenter")' - - julia --project ./docs/make.jl + - julia --project -e 'import Pkg; Pkg.add("Documenter"); Pkg.add("Coverage")' + - julia --project -e 'using Coverage; Codecov.submit(process_folder())' + - julia --project ./docs/make.jl notifications: email: true - slack: - rooms: - - secure: YydV4a8we+y26g/dlKsPFhXTOqJ7Q755YGOtGB19vo3JnG6QJ22DOQxF2+cMN3YMKRlwE1BhW6sB8Vd6BbfU4mPw8VXaSjnckLp4C48K6vxpX21gpeajFawaxAeyUPOp8KJZyKjNWGpz2UxNXv/f0V1iIGKs4ZJ0J3kg4/kMZvpOD+DiYoR6JN1xtu1yPbY67DogUB2LPfK1k/uiyaN166dzelVZtj1+bO3lKveBfB0KcI0HaUtRnnIzvZ65ZiuGJeTIFXoWl0MaKxlp4Nw5W3lz5pTbjmbBxtfpZBBz/70dNnUzbcpBzRiRRrkDu7bPT25ZgqS3ynHNaGH2LmLiv2T9/9A6Rqc6u/H8vfcBxNIbq+NHwUFQiNttLXtnsz20GOpVizfLvdM2s0/iZfP90ONtn/4o0/XGo5+dzYLMJvhh3osjKEx/O+fdNUVQJYHdg+aXnvEpcmwH7HGuDFJEIztwbjZNOoVgiDA3vf+H1bnFlKQV8mqmLygQ7k13iceQ2emEqxgaaaM52w4V3WiDDtdllILn6k3LVixERJv54OBzDjN0nvCXTYAFDnkTX8rvrZ+ieMb9CORd5ONBXl6MHWfbjFijdpQ973t2u7kxeOeggsk0S2vjSrw5IuuB+CfeyuLnKLurZ1H0fXywpUuhD22bhoSR0JZCJU3dCrjtmFc= - on_success: always - on_failure: always - template: - - "Repo `%{repository_slug}` *%{result}* build (<%{build_url}|#%{build_number}>) for commit (<%{compare_url}|%{commit}>) on branch `%{branch}`." - - "Execution time: *%{duration}*" - - "Message: %{message}" diff --git a/Project.toml b/Project.toml index 8e305c660..f36de4291 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "PorousMaterials" uuid = "68953c7c-a3c7-538e-83d3-73516288599e" authors = ["SimonEnsemble "] -version = "0.3.1" +version = "0.4.0" [deps] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" +FIGlet = "3064a664-84fe-4d92-92c7-ed492f3d8fae" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" @@ -18,13 +19,17 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Xtals = "ede5f01d-793e-4c47-9885-c447d1f18d6d" [compat] +FIGlet = "^0.2.1" +JLD2 = "^0.4.9" CSV = "^0.8.0" DataFrames = "^0.22.1" FileIO = "1.2.0" @@ -36,7 +41,8 @@ ProgressMeter = "1.2.0" Roots = "0.8.4" SpecialFunctions = "0.10.0" StatsBase = "0.32.0,0.33.8" -julia = "1.5,1.6" +Xtals = "0.3.0" +julia = "1.6" [extras] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" diff --git a/run_gcmc_long.sh b/run_gcmc_long.sh new file mode 100755 index 000000000..8b5d49566 --- /dev/null +++ b/run_gcmc_long.sh @@ -0,0 +1,15 @@ +#!/bin/env bash + +julia -p 24 -e ' + using Distributed + @everywhere begin + import Pkg + Pkg.activate(".") + end + #Pkg.test("PorousMaterials") + @info "Running test/gcmc_long.jl on $(length(workers())) workers." + @warn "This will take HOURS!" + cd("test") + include("gcmc_long.jl") + @info "Long GCMC tests complete!" +' diff --git a/src/PorousMaterials.jl b/src/PorousMaterials.jl index 52dcfe989..998d9875b 100644 --- a/src/PorousMaterials.jl +++ b/src/PorousMaterials.jl @@ -1,120 +1,19 @@ module PorousMaterials -using CSV -using DataFrames -using Roots # for fzero -using OffsetArrays # used for Ewald sum -using SpecialFunctions # for erfc -using StatsBase -using ProgressMeter -using Polynomials -using JLD2 -using Statistics -using Printf -using LinearAlgebra -using LightGraphs -using Distributed -using Optim -using PyCall - # import Base.push! - # - -# atoms are considered to overlap if this close. -const R²_OVERLAP = 0.1 # Units: Angstrom² - -""" - print_file_paths() - -print off paths where PorousMaterials.jl looks for input files and writes output files. -""" -function print_file_paths() - println("general data folder: ", PATH_TO_DATA) - println("\tcrystal structures (.cif, .cssr): ", PATH_TO_CRYSTALS) - println("\tforce field files (.csv): ", PATH_TO_FORCEFIELDS) - println("\tmolecule input files: ", PATH_TO_MOLECULES) - println("\tsimulation output files: ", PATH_TO_SIMS) - println("\tgrids (.cube): ", PATH_TO_GRIDS) -end +using Roots, OffsetArrays, SpecialFunctions, StatsBase, ProgressMeter, Polynomials, +JLD2, Statistics, Distributed, Optim, Printf, DataFrames, LightGraphs, CSV, LinearAlgebra -""" - set_path_to_data(path; print_paths=true) - -Set the path variables: `PATH_TO_DATA`, `PATH_TO_CRYSTALS`, `PATH_TO_FORCEFIELDS`, `PATH_TO_MOLECULES`, -`PATH_TO_GRIDS`, and `PATH_TO_SIMS`. The latter five paths are set relative to the root data path. - -# Arguments -- `path::String`: the absolute path to the root of the data directory. -- `print_paths::Bool`: set false to suppress printing of path values. -""" -function set_path_to_data(path::String; print_paths::Bool=true) - global PATH_TO_DATA = path - global PATH_TO_CRYSTALS = joinpath(PATH_TO_DATA, "crystals") - global PATH_TO_FORCEFIELDS = joinpath(PATH_TO_DATA, "forcefields") - global PATH_TO_MOLECULES = joinpath(PATH_TO_DATA, "molecules") - global PATH_TO_GRIDS = joinpath(PATH_TO_DATA, "grids") - global PATH_TO_SIMS = joinpath(PATH_TO_DATA, "simulations") - - if print_paths - print_file_paths() - end -end +# extend Xtals +using Reexport +@reexport using Xtals +import Xtals.Cart, Xtals.Frac, Xtals.write_xyz -""" - set_default_file_paths(print_paths=true) - -sets the default paths for where input files and some output files are stored. -to see current set up, call [`print_file_paths`](@ref) -""" -function set_default_file_paths(;print_paths::Bool=true) - # this is the main directory where crystal structures, forcefields, and molecules data is stored - global PATH_TO_DATA = joinpath(pwd(), "data") - - global PATH_TO_CRYSTALS = joinpath(PATH_TO_DATA, "crystals") - global PATH_TO_FORCEFIELDS = joinpath(PATH_TO_DATA, "forcefields") - global PATH_TO_MOLECULES = joinpath(PATH_TO_DATA, "molecules") - global PATH_TO_GRIDS = joinpath(PATH_TO_DATA, "grids") - global PATH_TO_SIMS = joinpath(PATH_TO_DATA, "simulations") - - if print_paths - print_file_paths() - end -end +# physical constants +const UNIV_GAS_CONST = 8.3144598e-5 # m³-bar/(K-mol) +const K_TO_KJ_PER_MOL = 8.3144598e-3 # kJ/(mol-K) +const BOLTZMANN = 1.38064852e7 # Boltmann constant (Pa-m3/K --> Pa-A3/K) -# this runs everytime porousmaterials is loaded, so if the user changes directory -# then the path_to_data will change as well -function __init__() - set_default_file_paths(print_paths=false) -end - - # """ - # set_tutorial_mode() - # - # places porousmaterials in "tutorial mode". it changes the `path_to_data` variable to - # the directory where the porousmaterials test data is stored. it can be used to - # follow examples shown in the readme. it displays a warning so that the user knows - # they are no longer using their own data. - # """ - # function set_tutorial_mode() - # new_path = joinpath(dirname(pathof(porousmaterials)), "..", "test", "data") - # if ! isdir(new_path) - # @error @sprintf("directory for testing data %s does not exist.\nnot entering tutorial mode.\n", new_path) - # else - # global path_to_data = new_path - # global path_to_crystals = joinpath(path_to_data, "crystals") - # global path_to_forcefields = joinpath(path_to_data, "forcefields") - # global path_to_molecules = joinpath(path_to_data, "molecules") - # global path_to_grids = joinpath(path_to_data, "grids") - # @warn "porousmaterials is now in tutorial mode. you have access to the testing data to experiment with porousmaterials.\nto reset to default file paths, call `set_default_file_paths()`\n" - # end - # end - # -include("matter.jl") -include("box.jl") -include("distance.jl") -include("misc.jl") include("isotherm_fitting.jl") -include("crystal.jl") -include("bonds.jl") include("forcefield.jl") include("molecule.jl") include("energy_utilities.jl") @@ -125,48 +24,23 @@ include("grid.jl") include("eos.jl") include("henry.jl") include("gcmc.jl") - # include("generic_rotations.jl") include("energy_min.jl") +include("atomic_masses.jl") +function __init__() + rc[:paths][:forcefields] = joinpath(rc[:paths][:data], "forcefields") + rc[:paths][:molecules] = joinpath(rc[:paths][:data], "molecules") + rc[:paths][:grids] = joinpath(rc[:paths][:data], "grids") + rc[:paths][:sims] = joinpath(rc[:paths][:data], "simulations") + append_atomic_masses() + @debug "Environment variables" rc +end -export - # porousmaterials.jl - print_file_paths, set_path_to_data, - - # matter.jl - Coords, Frac, Cart, Atoms, Charges, wrap!, neutral, net_charge, translate_by!, origin, - - # box.jl - Box, replicate, unit_cube, write_vtk, inside, fractional_coords, cartesian_coords, - - # distance.jl - nearest_image!, distance, overlap, remove_duplicates, - - # misc.jl - read_xyz, read_cpk_colors, write_xyz, read_atomic_masses, +export # isotherm_fitting.jl fit_adsorption_isotherm, - # crystal.jl - Crystal, strip_numbers_from_atom_labels!, assign_charges, - chemical_formula, molecular_weight, crystal_density, write_cif, has_charges, - apply_symmetry_operations, - - # bonds.jl - infer_bonds!, write_bond_information, BondingRule, bond_sanity_check, remove_bonds!, - infer_geometry_based_bonds!, cordero_covalent_atomic_radii, - - # construct_box, - # replicate, read_atomic_masses, charged, write_cif, assign_charges, - # is_symmetry_equal, apply_symmetry_rules, assert_p1_symmetry, infer_bonds!, - # remove_bonds!, compare_bonds_in_framework, wrap_atoms_to_unit_cell!, - # write_bond_information, is_bonded, default_bondingrules, has_same_sets_of_atoms_and_charges, - # distance, bond_sanity_check, - # - # # FrameworkOperations.jl - # partition_framework, subtract_atoms, - # # molecule.jl Molecule, translate_by!, translate_to!, random_rotation!, random_rotation_matrix, ion, distortion, @@ -189,20 +63,17 @@ export # Grid.jl apply_periodic_boundary_condition!, Grid, write_cube, read_cube, energy_grid, compute_accessibility_grid, accessible, - required_n_pts, xf_to_id, id_to_xf, update_density!, find_energy_minimum, - # + required_n_pts, xf_to_id, id_to_xf, update_density!, find_energy_minimum, origin, + # EOS.jl calculate_properties, PengRobinsonFluid, VdWFluid, # gcmc.jl μVT_sim, adsorption_isotherm, stepwise_adsorption_isotherm, μVT_output_filename, GCMCstats, MarkovCounts, isotherm_sim_results_to_dataframe, - # + # henry.jl henry_coefficient, henry_result_savename, - # - # # generic_rotations.jl - # rotation_matrix # energy_min.jl find_energy_minimum, find_energy_minimum_gridsearch diff --git a/src/atomic_masses.jl b/src/atomic_masses.jl new file mode 100644 index 000000000..170e58111 --- /dev/null +++ b/src/atomic_masses.jl @@ -0,0 +1,24 @@ +# add atom types to rc[:atomic_masses] for contextual LJ potentials +function append_atomic_masses() + rc[:atomic_masses] = merge(rc[:atomic_masses], Dict( + :N_in_N2 => 14.0067, + :CH2 => 14.025, + :CH3 => 15.035, + :CH4 => 16.04, + :C_b => 12.0107, + :C_tol => 12.0107, + :C_ac => 12.0107, + :C_RCOO => 12.0107, + :C_sp2 => 12.0107, + :C_sp3 => 12.0107, + :C_CO2 => 12.0107, + :H_H2S => 1.00794, + :H_b => 1.00794, + :O_RCOO => 15.9994, + :O_CO2 => 15.9994, + :O_zeo => 15.9994, + :S_H2S => 32.065, + :Si_zeo => 28.0855, + :X => 1.0 + )) +end \ No newline at end of file diff --git a/src/bonds.jl b/src/bonds.jl deleted file mode 100644 index ddef1d2a3..000000000 --- a/src/bonds.jl +++ /dev/null @@ -1,436 +0,0 @@ -""" - bonding_rule = BondingRule(:Ca, :O, 0.4, 2.0) - bonding_rules = [BondingRule(:H, :*, 0.4, 1.2), - BondingRule(:*, :*, 0.4, 1.9)] - -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 --`min_dist`: The minimum distance between the atoms for bonding to occur --`max_dist`: The maximum distance between the atoms for bonding to occur -""" -struct BondingRule - species_i::Symbol - species_j::Symbol - min_dist::Float64 - max_dist::Float64 -end - -""" - default_bondingrules = default_bondingrules() - -Returns the default bonding rules. Using `append!` and/or `prepend!` to add to the default bonding rules: - -# Example -``` -bond_rules = default_bondingrules() -prepend!(bond_rules, BondingRule(:Cu, :*, 0.1, 2.6)) -``` - -# Returns --`default_bondingrules::Array{BondingRule, 1}`: The default bonding rules: `[BondingRule(:*, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)]` -""" -default_bondingrules() = [BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)] - -""" - are_atoms_bonded = is_bonded(crystal, i, j, bonding_rules=[BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)], - include_bonds_across_periodic_boundaries=true) - -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 - 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 - periodic boundary when calculating bonds - -# Returns --`are_atoms_bonded::Bool`: Whether atoms `i` and `j` are bonded according to `bonding_rules` - -""" -function is_bonded(crystal::Crystal, i::Int64, j::Int64, - bonding_rules::Array{BondingRule, 1}=[BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)]; - include_bonds_across_periodic_boundaries::Bool=true) - species_i = crystal.atoms.species[i] - species_j = crystal.atoms.species[j] - - r = distance(crystal.atoms, crystal.box, i, j, include_bonds_across_periodic_boundaries) - - # 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) - species_match = true - 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) - species_match = true - end - - if species_match - # determine if the atoms are close enough to bond - if br.min_dist < r && br.max_dist > r - return true - else - return false # found relevant bonding rule, don't apply others - end - end - end - return false # no bonding rule applied -end - -""" - remove_bonds!(crystal) - -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) - end -end - -""" - infer_bonds!(crystal, include_bonds_across_periodic_boundaries, - bonding_rules=[BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)]) - -Populate the bonds in the crystal object based on the bonding rules. If a -pair doesn't have a suitable rule then they will not be considered bonded. - -`:*` is considered a wildcard and can be substituted for any species. It is a -good idea to include a bonding rule between two `:*` to allow any atoms to bond -as long as they are close enough. - -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. -""" -function infer_bonds!(crystal::Crystal, include_bonds_across_periodic_boundaries::Bool, - bonding_rules::Array{BondingRule, 1}= - [BondingRule(:H, :*, 0.4, 1.2), BondingRule(:*, :*, 0.4, 1.9)]) - @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 - # loop over every unique pair of atoms - for j in i+1:crystal.atoms.n - if is_bonded(crystal, i, j, bonding_rules; include_bonds_across_periodic_boundaries=include_bonds_across_periodic_boundaries) - add_edge!(crystal.bonds, i, j) - end - end - end -end - -""" - atom_to_radius = cordero_covalent_atomic_radii() - -Create a dictionary with the Cordero covalent radius for each element - -# Returns --`atom_to_radius::Dict{Symbol, Float64}`: A dictionary with elements as keys and their respective cordero covalent radii as the values. -""" -function cordero_covalent_atomic_radii() - df = CSV.read(joinpath(PATH_TO_DATA, "covalent_radii.csv"), DataFrame, comment="#") - atom_to_radius = Dict{Symbol, Float64}() - for atom in eachrow(df) - atom_to_radius[Symbol(atom[:atom])] = atom[:covalent_radius_A] - end - return atom_to_radius -end - -""" - a = adjacency_matrix(crystal, apply_pbc) - -Compute the adjacency matrix `a` of the crystal, where `a[i, j]` is the -distance between atom `i` and `j`. This matrix is symmetric. If `apply_pbc = true`, -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 --`a::Array{Float64, 2}`: symmetric, square adjacency matrix with zeros on the diagonal -""" -function adjacency_matrix(crystal::Crystal, apply_pbc::Bool) - A = zeros(crystal.atoms.n, crystal.atoms.n) - for i = 1:crystal.atoms.n - for j = (i+1):crystal.atoms.n - A[i, j] = distance(crystal.atoms, crystal.box, i, j, apply_pbc) - A[j, i] = A[i, j] # symmetric - end - end - return A -end - -""" - ids_neighbors, xs, rs = neighborhood(crystal, i, r, am) - -Find and characterize the neighborhood of atom `i` in the crystal `crystal`. -A neighborhood is defined as all atoms within a distance `r` from atom `i`. -The adjacency matrix `am` is used to find the distances of all other atoms in the crystal from atom `i`. - -# Arguments --`crystal::Crystal`: crystal structure --`i::Int`: Index of the atom (in `crystal`) which the neighborhood is to be characterized. --`r::Float64`: The maximum distance the neighborhood will be characterized. --`am::Array{Float64, 2}`: The adjacency matrix, see [`adjacency_matrix`](@ref) - -# Returns --`ids_neighbors::Array{Int, 1}`: indices of `crystal.atoms` within the neighborhood of atom `i`. --`xs::Array{Array{Float64, 1}, 1}`: array of Cartesian positions of the atoms surrounding atom `i`. - The nearest image convention has been applied to find the nearest periodic image. Also, the coordinates of atom `i` - have been subtracted off from these coordinates so that atom `i` lies at the origin of this new coordinate system. - The first vector in `xs` is `[0, 0, 0]` corresponding to atom `i`. - The choice of type is for the Voronoi decomposition in Scipy. --`rs::Array{Float64, 1}`: list of distances of the neighboring atoms from atom `i`. -""" -function neighborhood(crystal::Crystal, i::Int, r::Float64, am::Array{Float64, 2}) - # get indices of atoms within a distance r of atom i - # the greater than zero part is to not include itself - ids_neighbors = findall((am[:, i] .> 0.0) .& (am[:, i] .< r)) - - # rs is the list of distance of these neighbors from atom i - rs = [am[i, id_n] for id_n in ids_neighbors] - @assert all(rs .< r) - - # xs is a list of Cartesian coords of the neighborhood - # coords of atom i are subtracted off - # first entry is coords of atom i, the center, the zero vector - # remaining entries are neighbors - # this list is useful to pass to Voronoi for getting Voronoi faces - # of the neighborhood. - xs = [[0.0, 0.0, 0.0]] # this way atom zero is itself - for j in ids_neighbors - # subtract off atom i, apply nearest image - xf = crystal.atoms.coords.xf[:, j] - crystal.atoms.coords.xf[:, i] - nearest_image!(xf) - x = crystal.box.f_to_c * xf - push!(xs, x) - end - return ids_neighbors, xs , rs -end - -""" - ids_shared_voro_face = _shared_voronoi_faces(ids_neighbors, xs) - -Of the neighboring atoms, find those that share a Voronoi face. - -# Arguments --`ids_neighbors::Array{Int, 1}`: indices of atoms within the neighborhood of a specific atom. --`xs::Array{Array{Float64, 1}, 1}`: array of Cartesian position of the atoms within the neighborhood of a specific atom, relative to the specific atom. - -# Returns --`ids_shared_voro_face::Array{Int, 1}`: indices of atoms that share a Voronoi face with a specific atom -""" -function _shared_voronoi_faces(ids_neighbors::Array{Int, 1}, xs::Array{Array{Float64, 1}, 1}) - scipy = pyimport("scipy.spatial") - # first element of xs is the point itself, the origin - @assert length(ids_neighbors) == (length(xs) - 1) - - voro = scipy.Voronoi(xs) - rps = voro.ridge_points # connection with atom zero are connection with atom i - ids_shared_voro_face = Int[] # corresponds to xs, not to atoms of crystal - for k = 1:size(rps)[1] - if sort(rps[k, :])[1] == 0 # a shared face with atom i! - push!(ids_shared_voro_face, sort(rps[k, :])[2]) - end - end - # zero based indexing in Scipy accounted for since xs[0] is origin, atom i. - return ids_neighbors[ids_shared_voro_face] -end - -""" - ids_bonded = bonded_atoms(crystal, i, am; r=6.0, tol=0.25) - -Returns the ids of atoms that are bonded to atom `i` by determining bonds using a Voronoi method - -# Arguments --`crystal::Crystal`: Crystal structure in which the bonded atoms will be determined --`i::Int`: Index of the atom we want to determine the bonds of --`am::Array{Float64, 2}`: The adjacency matrix, see [`adjacency_matrix`](@ref) --`r::Float64`: The maximum distance used to determine the neighborhood of atom `i` --`tol::Float64`: A tolerance used when determining if two atoms are connected - -# Returns --`ids_bonded::Array{Int, 1}`: A list of indices of atoms bonded to atom `i` -""" -function bonded_atoms(crystal::Crystal, i::Int, am::Array{Float64, 2}; r::Float64=6.0, - tol::Float64=0.25, covalent_radii::Union{Nothing, Dict{Symbol, Float64}}=nothing) - if isnothing(covalent_radii) - covalent_radii = cordero_covalent_atomic_radii() - end - species_i = crystal.atoms.species[i] - - ids_neighbors, xs, rs = neighborhood(crystal, i, r, am) - - ids_shared_voro_faces = _shared_voronoi_faces(ids_neighbors, xs) - - ids_bonded = Int[] - for j in ids_shared_voro_faces - species_j = crystal.atoms.species[j] - # sum of covalent radii - radii_sum = covalent_radii[species_j] + covalent_radii[species_i] - if am[i, j] < radii_sum + tol - push!(ids_bonded, j) - end - end - return ids_bonded -end - -""" - infer_geometry_based_bonds!(crystal, include_bonds_across_periodic_boundaries::Bool) - -Infers bonds by first finding which atoms share a Voronoi face, and then bond the atoms if the distance - between them is less than the sum of the covalent radius of the two atoms (plus a tolerance). - -# Arguments --`crystal::Crystal`: The crystal structure --`include_bonds_across_periodic_boundaries::Bool`: Whether to check across the periodic boundaries -""" -function infer_geometry_based_bonds!(crystal::Crystal, include_bonds_across_periodic_boundaries::Bool; - r::Float64=6.0, tol::Float64=0.25, - covalent_radii::Union{Nothing, Dict{Symbol, Float64}}=nothing) - @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) - am = adjacency_matrix(crystal, include_bonds_across_periodic_boundaries) - - for i = 1:crystal.atoms.n - for j in bonded_atoms(crystal, i, am; r=r, tol=tol, covalent_radii=covalent_radii) - add_edge!(crystal.bonds, i, j) - end - 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 - -if sanity checks fail, refer to [`write_bond_information`](@ref) to write a .vtk to visualize the bonds. - -Print warnings when sanity checks fail. -Return `true` if sanity checks pass, `false` otherwise. -""" -function bond_sanity_check(crystal::Crystal) - for a = 1:crystal.atoms.n - ns = neighbors(crystal.bonds, a) - # is the graph fully connected? - if length(ns) == 0 - @warn "atom $a = $(crystal.atoms.species[a]) in $(crystal.name) is not bonded to any other atom." - return false - end - # does hydrogen have only one bond? - if (crystal.atoms.species[a] == :H) && (length(ns) > 1) - @warn "hydrogen atom $a in $(crystal.name) is bonded to more than one atom!" - return false - end - # does carbon have greater than four bonds? - if (crystal.atoms.species[a] == :C) && (length(ns) > 4) - @warn "carbon atom $a in $(crystal.name) is bonded to more than four atoms!" - return false - end - end - return true -end - -# TODO remove? why is this needed? -""" - bonds_equal = compare_bonds_in_crystal(crystal1, crystal2, atol=0.0) - -Returns whether the bonds defined in crystal1 are the same as the bonds -defined in crystal2. It checks whether the atoms in the same positions -have the same bonds. - -# Arguments --`crystal1::Crystal`: The first crystal --`crystal2::Crystal`: The second crystal --`atol::Float64`: absolute tolerance for the comparison of coordinates in the crystal - -# Returns --`bonds_equal::Bool`: Wether the bonds in crystal1 and crystal2 are equal -""" -function compare_bonds_in_crystal(fi::Crystal, fj::Crystal; atol::Float64=0.0) - if ne(fi.bonds) != ne(fj.bonds) - return false - end - - num_in_common = 0 - for edge_i in collect(edges(fi.bonds)) - for edge_j in collect(edges(fj.bonds)) - # either the bond matches going src-src dst-dst - if (fi.atoms.species[edge_i.src] == fj.atoms.species[edge_j.src] && - fi.atoms.species[edge_i.dst] == fj.atoms.species[edge_j.dst] && - isapprox(fi.atoms.xf[:, edge_i.src], fj.atoms.xf[:, edge_j.src]; atol=atol) && - isapprox(fi.atoms.xf[:, edge_i.dst], fj.atoms.xf[:, edge_j.dst]; atol=atol)) || - # or the bond matches going src-dst dst-src - (fi.atoms.species[edge_i.src] == fj.atoms.species[edge_j.dst] && - fi.atoms.species[edge_i.dst] == fj.atoms.species[edge_j.src] && - isapprox(fi.atoms.xf[:, edge_i.src], fj.atoms.xf[:, edge_j.dst]; atol=atol) && - isapprox(fi.atoms.xf[:, edge_i.dst], fj.atoms.xf[:, edge_j.src]; atol=atol)) - num_in_common += 1 - break - end - end - end - return num_in_common == ne(fi.bonds) && num_in_common == ne(fj.bonds) -end - -""" - write_bond_information(crystal, filename) - write_bond_information(crystal, center_at_origin=false) - -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`: center the coordinates at the origin of the crystal -""" -function write_bond_information(crystal::Crystal, filename::String; center_at_origin::Bool=false) - 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) - end - if ! occursin(".vtk", filename) - filename *= ".vtk" - end - - vtk_file = open(filename, "w") - - @printf(vtk_file, "# vtk DataFile Version 2.0\n%s bond information\nASCII\nDATASET POLYDATA\nPOINTS %d double\n", crystal.name, nv(crystal.bonds)) - - for i = 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]))...) - else - @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", ne(crystal.bonds), 3 * ne(crystal.bonds)) - for edge in collect(edges(crystal.bonds)) - @printf(vtk_file, "2\t%d\t%d\n", edge.src - 1, edge.dst - 1) - end - close(vtk_file) - @printf("Saving bond information for crystal %s to %s.\n", crystal.name, joinpath(pwd(), filename)) -end - -write_bond_information(crystal::Crystal; center_at_origin::Bool=false) = write_bond_information(crystal, split(crystal.name, ".")[1] * "_bonds.vtk", center_at_origin=center_at_origin) - -# TODO remove bonds with atom i? diff --git a/src/box.jl b/src/box.jl deleted file mode 100644 index 00ab8034f..000000000 --- a/src/box.jl +++ /dev/null @@ -1,249 +0,0 @@ -@doc raw""" - box = Box(a, b, c, α, β, γ, volume, f_to_c, c_to_f, reciprocal_lattice) - box = Box(a, b, c, α, β, γ) - box = Box(a, b, c) # α=β=γ=π/2 assumed. - box = Box(f_to_c) - -Data structure to describe a unit cell box (Bravais lattice) and convert between -fractional and Cartesian coordinates. - -# Attributes -- `a,b,c::Float64`: unit cell dimensions (units: Angstroms) -- `α,β,γ::Float64`: unit cell angles (units: radians) -- `Ω::Float64`: volume of the unit cell (units: cubic Angtroms) -- `f_to_c::Array{Float64,2}`: the 3x3 transformation matrix used to map fractional -coordinates to cartesian coordinates. The columns of this matrix define the unit cell -axes. Columns are the vectors defining the unit cell box. units: Angstrom -- `c_to_f::Array{Float64,2}`: the 3x3 transformation matrix used to map Cartesian -coordinates to fractional coordinates. units: inverse Angstrom -- `reciprocal_lattice::Array{Float64, 2}`: the *rows* are the reciprocal lattice vectors. -This choice was made (instead of columns) for speed of Ewald Sums. -""" -struct Box - a::Float64 - b::Float64 - c::Float64 - - α::Float64 - β::Float64 - γ::Float64 - - Ω::Float64 - - f_to_c::Array{Float64, 2} - c_to_f::Array{Float64, 2} - - reciprocal_lattice::Array{Float64, 2} -end - -Base.broadcastable(b::Box) = Ref(b) - -# Given the fractional to Cartesian transformation matrix, compute the reciprocal lattice -# vectors, which are the rows of the returned matrix `r`. -function reciprocal_lattice(f_to_c::Array{Float64, 2}) - # the unit cell vectors are the columns of f_to_c - a₁ = f_to_c[:, 1] - a₂ = f_to_c[:, 2] - a₃ = f_to_c[:, 3] - - r = zeros(Float64, 3, 3) - r[1, :] = 2 * π * cross(a₂, a₃) / dot(a₁, cross(a₂, a₃)) - r[2, :] = 2 * π * cross(a₃, a₁) / dot(a₂, cross(a₃, a₁)) - r[3, :] = 2 * π * cross(a₁, a₂) / dot(a₃, cross(a₁, a₂)) - return r -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) - # unit cell volume (A³) - Ω = 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(γ) / Ω]] - # the columns of f_to_c are the unit cell axes - r = reciprocal_lattice(f_to_c) - - @assert f_to_c * c_to_f ≈ Matrix{Float64}(I, 3, 3) - @assert isapprox(r, 2.0 * π * inv(f_to_c)) - - return Box(a, b, c, α, β, γ, Ω, f_to_c, c_to_f, r) -end - -# Constructs a `Box` from the fractional to Cartesian coordinate transformation matrix. -# The columns of this matrix are the unit cell axes. Units must be in Å. -function Box(f_to_c::Array{Float64, 2}) - # unit cell volume (A³) is the determinant of the matrix - Ω = det(f_to_c) - - # unit cell dimensions are lengths of the unit cell axes (Å) - a = norm(f_to_c[:, 1]) - b = norm(f_to_c[:, 2]) - c = norm(f_to_c[:, 3]) - - # c_to_f is the inverse - c_to_f = inv(f_to_c) - - # angles (radians) - α = acos(dot(f_to_c[:, 2], f_to_c[:, 3]) / (b * c)) - β = acos(dot(f_to_c[:, 1], f_to_c[:, 3]) / (a * c)) - γ = acos(dot(f_to_c[:, 1], f_to_c[:, 2]) / (a * b)) - - # the columns of f_to_c are the unit cell axes - r = reciprocal_lattice(f_to_c) - - return Box(a, b, c, α, β, γ, Ω, f_to_c, c_to_f, r) -end - -""" - uc = unit_cube() - -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 - -""" - f_coords = Frac(c_coords, box) - atoms_f = Frac(atoms_c, box) - charges_f = Frac(charges_c, box) - -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)) -Frac(charges::Charges{Cart}, box::Box) = Charges(charges.q, Frac(charges.coords, box)) - -""" - c_coords = Cart(f_coords, box) - 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. -""" -Cart(coords::Frac, box::Box) = Cart(box.f_to_c * coords.xf) -Cart(atoms::Atoms{Frac}, box::Box) = Atoms(atoms.species, Cart(atoms.coords, box)) -Cart(charges::Charges{Frac}, box::Box) = Charges(charges.q, Cart(charges.coords, box)) - -""" - new_box = replicate(original_box, repfactors) - -Replicates a `Box` in positive directions to construct a new `Box` representing a supercell. -The `original_box` is replicated according to the factors in `repfactors`. -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 - -# Returns -- `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.γ) -end - -""" - write_vtk(box, filename; verbose=true, center_at_origin=false) - write_vtk(framework) - -Write a `Box` to a .vtk file for visualizing e.g. the unit cell boundary of a crystal. -If a `Framework` is passed, the `Box` of that framework is written to a file that is the -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) -- `framework::Framework`: A framework containing the crystal structure information -- `center_at_origin::Bool`: center box at origin if true. if false, the origin is the corner of the box. -- `verbose::Bool`: Optional argument. If true (default), the output filename is printed to the console. -""" -function write_vtk(box::Box, filename::AbstractString; verbose::Bool=true, - center_at_origin::Bool=false) - if ! occursin(".vtk", filename) - filename *= ".vtk" - end - vtk_file = open(filename, "w") - - @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 - x_shift = box.f_to_c * [0.5, 0.5, 0.5] - end - - # write points on boundary of unit cell - for i = 0:1 - for j = 0:1 - for k = 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]) - 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") - close(vtk_file) - if verbose - println("See ", filename) - end -end - -""" - inside_box = inside(frac_coords) # true or false - inside_box = inside(cart_coords, box) - inside_box = inside(molecule, box) # in Cartesian - inside_box = inside(molecule) # in fractional - inside_box = inside(crystal) # in fractional - -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) -""" -inside(coords::Frac) = all(coords.xf .<= 1.0) && all(coords.xf .>= 0.0) -inside(coords::Cart, box::Box) = inside(Frac(coords, box)) - -# documented in matter.jl -translate_by!(coords::Frac, dx::Cart, box::Box) = translate_by!(coords, Frac(dx, box)) -translate_by!(coords::Cart, dx::Frac, box::Box) = translate_by!(coords, Cart(dx, box)) - -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, "\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)) -end diff --git a/src/crystal.jl b/src/crystal.jl deleted file mode 100644 index 86b1929b9..000000000 --- a/src/crystal.jl +++ /dev/null @@ -1,1077 +0,0 @@ -""" - SymmetryInfo(symmetry, space_group, is_p1) - -# Attributes -- `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 - 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. -""" -struct SymmetryInfo - operations::Array{String, 2} - space_group::String - is_p1::Bool -end -SymmetryInfo() = SymmetryInfo([Array{String, 2}(undef, 3, 0) ["x", "y", "z"]], "P1", true) # default - -struct Crystal - name::String - box::Box - atoms::Atoms{Frac} - charges::Charges{Frac} - bonds::SimpleGraph - symmetry::SymmetryInfo -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, SimpleGraph(atoms.n), SymmetryInfo()) - -const NET_CHARGE_TOL = 1e-4 # net charge tolerance - -""" - crystal = Crystal(filename; - check_neutrality=true, net_charge_tol=1e-4, - check_overlap=true, overlap_tol=0.1, - convert_to_p1=true, read_bonds_from_file=false, wrap_coords=true, - include_zero_charges=false, - remove_duplicates=false, - species_col=["_atom_site_label", "_atom_site_type_symbol"] - ) # read from file - - crystal = Crystal(name, box, atoms, charges) # construct from matter, no bonds, P1-symmetry assumed - -Read a crystal structure file (.cif or .cssr) and populate a `Crystal` data structure, -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 - 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 - 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. - 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. - -# Returns -- `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::SimpleGraph`: Unweighted, undirected graph showing all of the atoms - that are bonded within the crystal -- `symmetry::SymmetryInfo`: symmetry inforomation -- `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. -""" -function Crystal(filename::String; - check_neutrality::Bool=true, - net_charge_tol::Float64=NET_CHARGE_TOL, - check_overlap::Bool=true, - overlap_tol::Float64=0.1, - convert_to_p1::Bool=true, - read_bonds_from_file::Bool=false, - wrap_coords::Bool=true, - include_zero_charges::Bool=false, - remove_duplicates::Bool=false, - species_col::Array{String, 1}=["_atom_site_label", "_atom_site_type_symbol"]) - # Read file extension. Ensure we can read the file type - extension = split(filename, ".")[end] - if ! (extension in ["cif", "cssr"]) - error("I can only read .cif or .cssr crystal structure files.") - end - - # read all lines of crystal structure file - _f = open(joinpath(PATH_TO_CRYSTALS, filename), "r") - lines = readlines(_f) - close(_f) - - # Initialize arrays. We'll populate them when reading through the crystal structure file. - charge_values = Array{Float64, 1}() - species = Array{Symbol, 1}() - xf = Array{Float64, 1}() - yf = Array{Float64, 1}() - zf = Array{Float64, 1}() - coords = Array{Float64, 2}(undef, 3, 0) - # default for symmetry rules is P1. - # These will be overwritten if the user chooses to read in non-P1 - operations = Array{String, 2}(undef, 3, 0) - # creating empty SimpleGraph, might not have any information read in - bonds = SimpleGraph() - # used for remembering whether fractional/cartesian coordinates are read in - # placed here so it will be defined for the if-stmt after the box is defined - fractional = false - cartesian = false - # used for determining if the crystal is in P1 symmetry for simulations - is_p1 = false - space_group = "" - - - # Start of .cif reader - ################################### - # CIF READER - ################################### - if extension == "cif" - data = Dict{AbstractString, Float64}() - loop_starts = -1 - i = 1 - # used for reading in symmetry options and replications - symmetry_info = false - atom_info = false - label_num_to_idx = Dict{AbstractString, Int}() - fractional = false - cartesian = false - # find later. - species_column = "null" - - while i <= length(lines) - line = split(lines[i]) - # Skip empty lines - if length(line) == 0 - i += 1 - continue - end - - # Make sure the space group is P1 - if line[1] == "_symmetry_space_group_name_H-M" - # 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" - # simplify by only having one P1 space_group name - space_group = "P1" - is_p1 = true - end - end - - # checking for information about atom sites and symmetry - if line[1] == "loop_" - # creating dictionary of column names to determine what should be done - #atom_column_name = "" - # name_to_column is a dictionary that e.g. returns which column contains x fractional coord - # use example: name_to_column["_atom_site_fract_x"] gives 3 - name_to_column = Dict{AbstractString, Int}() - - i += 1 - loop_starts = i - while length(split(lines[i])) == 1 && split(lines[i])[1][1] == '_' - name_to_column[split(lines[i])[1]] = i + 1 - loop_starts - # iterate to next line in file - 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") - # 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 - - # 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 - found_species_col = false - for col in species_col - if col ∈ keys(name_to_column) - species_column = col - found_species_col = true - break - end - end - 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 - - # ===================== - # SYMMETRY READER - # ===================== - if haskey(name_to_column, "_symmetry_equiv_pos_as_xyz") - symmetry_info = true - - symmetry_count = 0 - # CSD stores symmetry as one column in a string that ends - # up getting split on the spaces between commas (i.e. its - # 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]) - line = lines[i] - sym_funcs = split(line, [' ', ',', ''', '"'], keepempty=false) - if length(sym_funcs) == 0 - break - end - symmetry_count += 1 - - # store as strings so it can be written out later - new_sym_rule = Array{AbstractString, 1}(undef, 3) - - sym_start = name_to_column["_symmetry_equiv_pos_as_xyz"] - 1 - for j = 1:3 - new_sym_rule[j] = sym_funcs[j + sym_start] - end - - operations = [operations new_sym_rule] - - i += 1 - end - - @assert symmetry_count == size(operations, 2) "number of symmetry rules must match the count" - - # finish reading in symmetry information, skip to next - # iteration of outer while-loop - continue - # ===================== - # FRACTIONAL READER - # ===================== - elseif fractional && ! atom_info - atom_info = true - - 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)]] - # 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"]])) - else - push!(charge_values, 0.0) - 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) - end - # iterate to next line in file - i += 1 - end - # set up graph of correct size - bonds = SimpleGraph(length(species)) - # finish reading in atom_site information, skip to next - # iteration of outer while-loop - # prevents skipping a line after finishing reading atoms - continue - # ===================== - # CARTESIAN READER - # ===================== - elseif cartesian && ! atom_info - atom_info = true - - 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])]] - # 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"]])) - else - push!(charge_values, 0.0) - 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) - end - # iterate to next line in file - i += 1 - end - # set up graph of correct size - bonds = SimpleGraph(length(species)) - # finish reading in atom_site information, skip to next - # iteration of outer while-loop - # prevents skipping a line after finishing reading atoms - continue - # ===================== - # 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) - 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) - - # iterate to next line in file - i += 1 - end - - # skip to next iteration in outer while loop - continue - end - end - - # 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]) - 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 - end - end - - i += 1 - end # End loop over lines - - if !atom_info - error("Could not find _atom_site* after loop_ in .cif file\n") - 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)) - end - - a = data["a"] - b = data["b"] - c = data["c"] - α = data["alpha"] - β = data["beta"] - γ = data["gamma"] - - # redo coordinates if they were read in cartesian - 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 - ################################### - elseif extension == "cssr" - # First line contains unit cell lenghts - line = split(lines[1]) - a = parse(Float64, line[1]) - b = parse(Float64, line[2]) - c = parse(Float64, line[3]) - - # Second line contains unit cell angles - line = split(lines[2]) - α = parse(Float64, line[1]) * pi / 180.0 - β = parse(Float64, line[2]) * pi / 180.0 - γ = parse(Float64, line[3]) * pi / 180.0 - - n_atoms = parse(Int, split(lines[3])[1]) - bonds = SimpleGraph(n_atoms) - - # Read in atoms and fractional coordinates - for i = 1:n_atoms - line = split(lines[4 + i]) - push!(species, Symbol(line[2])) - - push!(xf, parse(Float64, line[3])) - push!(yf, parse(Float64, line[4])) - push!(zf, parse(Float64, line[5])) - - push!(charge_values, parse(Float64, line[14])) - end - - for i = 1:n_atoms - coords = [ coords [xf[i], yf[i], zf[i]] ] - end - - # add P1 symmetry rules for consistency - operations = [operations ["x", "y", "z"]] - is_p1 = true - space_group = "P1" - end - - # Construct the unit cell box - box = Box(a, b, c, α, β, γ) - # construct atoms attribute of crystal - atoms = Atoms(species, Frac(coords)) - # construct charges attribute of crystal - if ! include_zero_charges - # include only nonzero charges - idx_nz = charge_values .!= 0.0 - charges = Charges(charge_values[idx_nz], Frac(coords[:, idx_nz])) - else - # include all charges, even if some are zero. - charges = Charges(charge_values, Frac(coords)) - end - - 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 afrain from this, pass `convert_to_p1=false` to the `Crystal` constructor.\n", - crystal.name, crystal.symmetry.space_group) - crystal = apply_symmetry_operations(crystal) - end - - if wrap_coords - wrap!(crystal) # do before checking overlap! - end - - if remove_duplicates - crystal = remove_duplicate_atoms_and_charges(crystal) - 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_charge_neutrality=false or increasing the - net charge tolerance `net_charge_tol`\n", crystal.name, net_charge(crystal))) - end - end - - if check_overlap - _check_overlap(crystal) - end - - return crystal -end - -overlap(crystal::Crystal) = overlap(crystal.atoms.coords, crystal.box, true) - -function _check_overlap(crystal::Crystal) - overlap_flag, overlapping_pairs = overlap(crystal) - if overlap_flag - for (i, j) in overlapping_pairs - @warn @sprintf("atom %d (%s) and %d (%s) are overlapping\n", i, crystal.atoms.species[i], - j, crystal.atoms.species[j]) - end - error(crystal.name * " has overlapping pairs of atoms! - (this occurs often when applying symmetry rules, if your structure was not in P1 symmetry to begin with.) - pass `check_overlap=false` then run `overlap(crystal)` to find pairs of overlapping atoms for inspection. - or pass `remove_duplicates=true` to the `Crystal` constructor to automatically remove the duplicate atoms and charges. - you should then visualize your .cif to make sure it is proper.`\n") - end -end - -# documented in matter.jl -function wrap!(crystal::Crystal) - wrap!(crystal.atoms.coords) - wrap!(crystal.charges.coords) -end - -# documented in matter.jl -net_charge(crystal::Crystal) = net_charge(crystal.charges) - -""" - replicated_crystal = replicate(crystal, repfactors) - -replicate the atoms and charges in a `Crystal` in positive directions to construct a new -`Crystal`. Note `replicate(crystal, (1, 1, 1))` returns the same `Crystal`. the fractional -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). - -# returns -- `replicated_frame::Crystal`: replicated crystal -""" -function replicate(crystal::Crystal, repfactors::Tuple{Int, Int, Int}) - if ne(crystal.bonds) != 0 - error("the crystal " * crystal.name * " has assigned bonds. to replicate, remove - its bonds with `remove_bonds!(crystal)`. then use `infer_bonds(crystal)` to - reassign the bonds") - end - - assert_P1_symmetry(crystal) - - n_atoms = crystal.atoms.n * prod(repfactors) - n_charges = crystal.charges.n * prod(repfactors) - - box = replicate(crystal.box, repfactors) - atoms = Atoms{Frac}(n_atoms) - charges = Charges{Frac}(n_charges) - - atom_counter = 0 - charge_counter = 0 - for ra = 0:(repfactors[1] - 1), rb = 0:(repfactors[2] - 1), rc = 0:(repfactors[3] - 1) - xf_shift = 1.0 * [ra, rb, rc] - - # replicate atoms - for i = 1:crystal.atoms.n - atom_counter += 1 - - atoms.species[atom_counter] = crystal.atoms.species[i] - - xf = crystal.atoms.coords.xf[:, i] + xf_shift - atoms.coords.xf[:, atom_counter] = xf ./ repfactors - end - - # replicate charges - for i = 1:crystal.charges.n - charge_counter += 1 - - charges.q[charge_counter] = crystal.charges.q[i] - - xf = crystal.charges.coords.xf[:, i] + xf_shift - charges.coords.xf[:, charge_counter] = xf ./ repfactors - end - end - - return Crystal(crystal.name, box, atoms, charges, SimpleGraph(n_atoms), crystal.symmetry) -end - -# doc string in Misc.jl -xyz_filename(crystal::Crystal) = replace(replace(crystal.name, ".cif" => ""), ".cssr" => "") * ".xyz" -function write_xyz(crystal::Crystal; comment::AbstractString="", center_at_origin::Bool=false) - filename = xyz_filename(crystal) - 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 - write_xyz(atoms, filename, comment=comment) - else - write_xyz(atoms, filename, comment=comment) - end -end - -# 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 - -`true` if any only if `crystal::Crystal`/`molecule::Molecule` has point charges. -""" -has_charges(crystal::Crystal) = crystal.charges.n > 0 - -""" - strip_numbers_from_atom_labels!(crystal) - -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 - -# Arguments -- `crystal::Crystal`: The crystal containing the crystal structure information -""" -function strip_numbers_from_atom_labels!(crystal::Crystal) - for i = 1:crystal.atoms.n - # atom species in string format - species = string(crystal.atoms.species[i]) - for j = 1:length(species) - if ! isletter(species[j]) - crystal.atoms.species[i] = Symbol(species[1:j-1]) - break - end - end - end -end - -vtk_filename(crystal::Crystal) = replace(replace(crystal.name, ".cif" => ""), ".cssr" => "") * ".vtk" -write_vtk(crystal::Crystal; center_at_origin::Bool=false) = write_vtk(crystal.box, vtk_filename(crystal), center_at_origin=center_at_origin) - -""" - formula = chemical_formula(crystal, verbose=false) - -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 - -# Returns -- `formula::Dict{Symbol, Int}`: A dictionary with the irreducible chemical formula of a crystal structure -""" -function chemical_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 - atom_counts[crystal.atoms.species[i]] += 1 - end - - # get greatest common divisor - gcd_ = gcd([k for k in values(atom_counts)]) - - # turn into irreducible chemical formula - for atom in keys(atom_counts) - atom_counts[atom] = atom_counts[atom] / gcd_ - end - - # print result - if verbose - @printf("Chemical formula of %s:\n\t", crystal.name) - for (atom, formula_unit) in atom_counts - @printf("%s_%d", string(atom), formula_unit) - end - @printf("\n") - end - - return atom_counts -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 - -# Returns -- `mass_of_crystal::Float64`: The molecular weight of a unit cell of the crystal in amu -""" -function molecular_weight(crystal::Crystal) - atomic_masses = read_atomic_masses() - - mass = 0.0 - for i = 1:crystal.atoms.n - mass += atomic_masses[crystal.atoms.species[i]] - end - - return mass # amu -end - -""" - ρ = crystal_density(crystal) # kg/m³ - -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 - -# Returns -- `ρ::Float64`: The crystal density of a crystal in kg/m³ -""" -crystal_density(crystal::Crystal) = molecular_weight(crystal) / crystal.box.Ω * 1660.53892 # --> kg/m3 - -""" - simulation_ready_crystal = apply_symmetry_operations(non_p1_crystal) - -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 - -# Returns -- `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) - if crystal.symmetry.is_p1 - return crystal - end - - nb_symmetry_ops = size(crystal.symmetry.operations, 2) - - n_atoms = crystal.atoms.n * nb_symmetry_ops - atoms = Atoms{Frac}(n_atoms) - - n_charges = crystal.charges.n * nb_symmetry_ops - charges = Charges{Frac}(n_charges) - - # for each symmetry rule - 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 - # 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] - end - - # loop over all charges in lower level symmetry - 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[c] = crystal.charges.q[c] - 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 - - # """ - # symmetry_equal = is_symmetry_equal(crystal1.symmetry, crystal2.symmetry) - # - # Returns true if both symmetry rules can create the same set from the same set - # of coordinates. Returns false if they don't contain the same number of rules or - # if they create different sets of points. - # - # # Arguments - # - `sym1::Array{AbstractString, 2}`: Array of strings that represent - # symmetry operations - # - `sym2::Array{AbstractString, 2}`: Array of strings that represent - # symmetry operations - # - # # Returns - # - `is_equal::Bool`: True if they are the same set of symmetry rules - # False if they are different - # """ - # function is_symmetry_equal(sym1::Array{AbstractString, 2}, sym2::Array{AbstractString, 2}) - # # need same number of symmetry operations - # if size(sym1, 2) != size(sym2, 2) - # return false - # end - # # define a test array that operations will be performed on - # test_array = [0.0 0.25 0.0 0.0 0.0 0.25 0.25 0.25; - # 0.0 0.0 0.25 0.0 0.25 0.0 0.25 0.25; - # 0.0 0.0 0.0 0.25 0.25 0.25 0.25 0.25] - # # set up both arrays for storing replicated coords - # sym1_applied_to_test = Array{Float64, 2}(undef, 3, 0) - # sym2_applied_to_test = Array{Float64, 2}(undef, 3, 0) - # - # # loop over all positions in the test_array - # for i in 1:size(test_array, 2) - # # loop over f1 symmetry rules - # for j in 1:size(sym1, 2) - # sym1_applied_to_test = [sym1_applied_to_test [Base.invokelatest.( - # eval(Meta.parse("(x, y, z) -> " * sym1[k, j])), test_array[:, i]...) for k in 1:3]] - # end - # # loop over f2 symmetry rules - # for j in 1:size(sym2, 2) - # sym2_applied_to_test = [sym2_applied_to_test [Base.invokelatest.( - # eval(Meta.parse("(x, y, z) -> " * sym2[k, j])), test_array[:, i]...) for k in 1:3]] - # end - # end - # - # # convert to sets for using issetequal, symmetry rules might be in a a different order - # sym1_set = Set([sym1_applied_to_test[:, i] for i in 1:size(sym1_applied_to_test, 2)]) - # sym2_set = Set([sym2_applied_to_test[:, i] for i in 1:size(sym2_applied_to_test, 2)]) - # - # # return if the sets of coords are equal - # return issetequal(sym1_set, sym2_set) - # 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 - 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 a `crystal::Crystal` to a .cif file with `filename::AbstractString`. If `filename` does -not include the .cif extension, it will automatically be added. the `fractional_coords` flag -allows us to write either fractional or Cartesian coordinates. -""" -function write_cif(crystal::Crystal, filename::AbstractString; 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) - error("write_cif needs coords of atoms and charges to correspond.") - end - end - - # TODO is this labeling necessary for the bonds, arthur? - # create dictionary for tracking label numbers - label_numbers = Dict{Symbol, Int}() - for atom in crystal.atoms.species - 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) - filename *= ".cif" - end - cif_file = open(filename, "w") - # first line should be data_xtalname_PM - if crystal.name == "" - @printf(cif_file, "data_PM\n") - else - # don't include file extension! - @printf(cif_file, "data_%s_PM\n", split(crystal.name, ".")[1]) - end - - @printf(cif_file, "_symmetry_space_group_name_H-M\t'%s'\n", crystal.symmetry.space_group) - - @printf(cif_file, "_cell_length_a\t%f\n", crystal.box.a) - @printf(cif_file, "_cell_length_b\t%f\n", crystal.box.b) - @printf(cif_file, "_cell_length_c\t%f\n", crystal.box.c) - - @printf(cif_file, "_cell_angle_alpha\t%f\n", crystal.box.α * 180.0 / pi) - @printf(cif_file, "_cell_angle_beta\t%f\n", crystal.box.β * 180.0 / pi) - @printf(cif_file, "_cell_angle_gamma\t%f\n", crystal.box.γ * 180.0 / pi) - - @printf(cif_file, "_symmetry_Int_Tables_number 1\n\n") - @printf(cif_file, "loop_\n_symmetry_equiv_pos_as_xyz\n") - for i in 1:size(crystal.symmetry.operations, 2) - @printf(cif_file, "'%s,%s,%s'\n", crystal.symmetry.operations[:, i]...) - end - @printf(cif_file, "\n") - - @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") - else - @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) - @printf(cif_file, "_atom_site_charge\n") - # if crystal will not be charge neutral to a 1e-5 tolerance when loading it - # into PorousMaterials.jl, then write higher-precision charges - if abs(sum(round.(crystal.charges.q, digits=6))) > NET_CHARGE_TOL - @info "writing high-precision charges for " * filename * ".\n" - high_precision_charges = true - end - end - - idx_to_label = Array{AbstractString, 1}(undef, crystal.atoms.n) - for i = 1:crystal.atoms.n - # print label and type symbol - @printf(cif_file, "%s\t%s\t", 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]]) - # increment label - label_numbers[crystal.atoms.species[i]] += 1 - if fractional_coords - @printf(cif_file, "%f\t%f\t%f", crystal.atoms.coords.xf[:, i]...) - else - @printf(cif_file, "%f\t%f\t%f", (crystal.box.f_to_c * crystal.atoms.coords.xf[:, i])...) - end - if has_charges(crystal) - if high_precision_charges - @printf(cif_file, "\t%.10f\n", crystal.charges.q[i]) - else - @printf(cif_file, "\t%f\n", crystal.charges.q[i]) - end - else - @printf(cif_file, "\n") - end - end - - # only print bond information if it is in the crystal - if ne(crystal.bonds) > 0 - 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") - - for edge in collect(edges(crystal.bonds)) - dxf = crystal.atoms.coords.xf[:, edge.src] - crystal.atoms.coords.xf[:, edge.dst] - nearest_image!(dxf) - @printf(cif_file, "%s\t%s\t%0.5f\n", idx_to_label[edge.src], idx_to_label[edge.dst], - norm(dxf)) - end - end - close(cif_file) -end - -write_cif(crystal::Crystal) = write_cif(crystal, crystal.name) - -""" - crystal = remove_duplicate_atoms_and_charges(crystal, r_tol=0.1, q_tol=0.0001, verbose=false) - -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 - -# returns -- `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) - 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, SimpleGraph(crystal.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) - -assign charges to the atoms present in the crystal based on atom type. -pass a dictionary `species_to_charge` that maps atomic species to a charge. - -if the crystal already has charges, the charges are removed and new charges are added. a warning is thrown if this is the case. - -checks for charge neutrality in the end. - -returns a new crystal. - -# Examples - -``` -species_to_charge = Dict(:Ca => 2.0, :C => 1.0, :H => -1.0) -crystal_with_charges = assign_charges(crystal, species_to_charge, 1e-7) -crystal_with_charges = assign_charges(crystal, species_to_charge) # tol 1e-5 default -``` - -# 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 -""" -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) - 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) - - # 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)) - end - - return new_crystal -end - -function Base.show(io::IO, crystal::Crystal) - println(io, "Name: ", crystal.name) - println(io, crystal.box) - @printf(io, "\t# atoms = %d\n", crystal.atoms.n) - @printf(io, "\t# charges = %d\n", crystal.charges.n) - println(io, "\tchemical formula: ", chemical_formula(crystal)) - @printf(io, "\tspace Group: %s\n", crystal.symmetry.space_group) - @printf(io, "\tsymmetry Operations:\n") - for i in 1:size(crystal.symmetry.operations, 2) - @printf(io, "\t\t'%s, %s, %s'\n", crystal.symmetry.operations[:, i]...) - end -end - -function Base.isapprox(c1::Crystal, c2::Crystal; atol::Real=0.0) - # name not included - box_flag = isapprox(c1.box, c2.box) - 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) - sym1 = c1.symmetry - sym2 = c2.symmetry - symmetry_flag = (sym1.is_p1 == sym2.is_p1) && (sym1.space_group == sym1.space_group) && (sym1.operations == sym2.operations) - return box_flag && charges_flag && atoms_flag && symmetry_flag -end - -function Base.getindex(crystal::Crystal, ids::Union{BitArray{1}, Array{Int, 1}, - UnitRange{Int}}) -# if ne(crystal.bonds) != 0 -# error("indexing a crystal with bonds is not supported.") -# end - - old_to_new = Dict(ids[i] => i for i in 1:length(ids)) - bonds = SimpleGraph(length(ids)) - for edge in collect(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]) - end - end - - if crystal.charges.n == 0 - 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) - else - error("for getindex(crystal), crystal must have 0 charges or an equal number of charges and atoms - that share coordinates") - end - - return crystal -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") - end -end - -function Base.:+(crystals::Crystal...; check_overlap::Bool=true) - crystal = deepcopy(crystals[1]) - for (i, f) in enumerate(crystals) - if i == 1 - continue - end - @assert isapprox(crystal.box, f.box) @sprintf("Crystal %s has a different box\n", f.name) - # @assert is_symmetry_equal(crystal.symmetry, f.symmetry) @sprintf("Crystal %s has different symmetry rules\n", f.name) - @assert crystal.symmetry != f.symmetry @sprintf("Crystal %s has different symmetry rules\n", f.name) - @assert crystal.symmetry.space_group == f.symmetry.space_group - - atoms = crystal.atoms + f.atoms - charges = crystal.charges + f.charges - - nf_n_atoms = crystal.atoms.n - add_vertices!(crystal.bonds, nf_n_atoms) - for edge in collect(edges(f.bonds)) - add_edge!(crystal.bonds, nf_n_atoms + edge.src, nf_n_atoms + edge.dst) - end - - crystal = Crystal(split(crystal.name, ".")[1] * "_" * split(f.name, ".")[1], - crystal.box, atoms, charges, crystal.bonds, - crystal.symmetry) - end - - if check_overlap - _check_overlap(crystal) - end - - return crystal -end diff --git a/src/distance.jl b/src/distance.jl deleted file mode 100644 index cacdfb9c0..000000000 --- a/src/distance.jl +++ /dev/null @@ -1,157 +0,0 @@ -""" - nearest_image!(dxf) - -Applies the nearest image convention on a vector `dxf` between two atoms in fractional -space; modifies `dxf` for nearest image convention. Fractional coordinates here fall in -[0, 1] so that the box is [0, 1]^3 in fractional space. - -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 -""" -@inline function 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) - r = distance(atoms, box, i, j, apply_pbc) # atoms i and j - r = distance(charges, box, i, j, apply_pbc) # atoms i and j - -calculate the (Cartesian) distance between particles `i` and `j`. - -apply periodic boundary conditions if and only if `apply_pbc` is `true`. - -# arguments -- `coords::Coords`: the coordinates (`Frac>:Coords` or `Cart>:Coords`) -- `atoms::Atoms`: atoms -- `charges::charges`: atoms -- `box::Box`: unit cell information -- `i::Int`: index of the first particle -- `j::Int`: Index of the second particle -- `apply_pbc::Bool`: `true` if we wish to apply periodic boundary conditions, `false` otherwise -""" -function distance(coords::Frac, box::Box, i::Int, j::Int, apply_pbc::Bool) - dxf = coords.xf[:, i] - coords.xf[:, j] - if apply_pbc - nearest_image!(dxf) - end - return norm(box.f_to_c * dxf) -end - -function distance(coords::Cart, box::Box, i::Int, j::Int, apply_pbc::Bool) - dx = coords.x[:, i] - coords.x[:, j] - if apply_pbc - dxf = box.c_to_f * dx - nearest_image!(dxf) - return norm(box.f_to_c * dxf) - else - return norm(dx) - end -end - -distance(atoms::Atoms, box::Box, i::Int, j::Int, apply_pbc::Bool) = distance(atoms.coords, box, i, j, apply_pbc) -distance(charges::Charges, box::Box, i::Int, j::Int, apply_pbc::Bool) = distance(charges.coords, box, i, j, apply_pbc) - -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 - pd[i, j] = distance(coords, box, i, j, apply_pbc) - if i > j - 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) - -""" - overlap_flag, overlap_pairs = overlap(frac_coords, box, apply_pbc; tol=0.1) - overlap_flag, overlap_pairs = overlap(crystal) - -determine if any coordinates overlap. here, two coordinates are defined to overlap if their (Cartesian) distance -is less than `tol`. - -# 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 - -# 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)]` -""" -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 - r = distance(coords, box, i, j, apply_pbc) - if r < tol - push!(overlap_ids, (i, j)) - overlap_flag = true - end - end - end - 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) - -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 - -# 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 -""" -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 - if isa(ac, Atoms{Frac}) - # if different species, not duplicates. - if ac.species[i] != ac.species[j] - continue - end - elseif isa(ac, Charges{Frac}) - # if different value of charge, not duplicates. - if ! isapprox(ac.q[i], ac.q[j], atol=q_tol) - continue - end - end - # are they overlapping? - r = distance(ac.coords, box, i, j, apply_pbc) - if r < r_tol - ids_keep[j] = false - end - end - end - return ac[ids_keep] -end diff --git a/src/energy_min.jl b/src/energy_min.jl index 3b191fbe4..5a9258961 100644 --- a/src/energy_min.jl +++ b/src/energy_min.jl @@ -66,19 +66,16 @@ perform an [`energy_grid`](@ref) calculation and, via a grid search, find the mi - `xtal::Crystal`: The crystal being investigated - `molecule::Molecule{Cart}`: The molecule used to probe energy surface - `ljff::LJForceField`: The force field used to calculate interaction energies -- `n_pts::Tuple{Int, Int, Int}=(50,50,50)`: Number of grid points in each fractional coordinate dimension, including endpoints (0, 1) +- `resolution::Union{Float64, Tuple{Int, Int, Int}}=1.0`: maximum distance between grid points, in Å, or a tuple specifying the number of grid points in each dimension. # Returns - `minimized_molecule::Molecule{Frac}`: the molecule at its minimum energy position - `min_energy::Float64`: the associated minimum molecule-crystal interaciton energy (kJ/mol) """ -function find_energy_minimum_gridsearch(xtal::Crystal, - molecule::Molecule{Cart}, - ljff::LJForceField; - n_pts::Tuple{Int, Int, Int}=(50, 50, 50) - ) +function find_energy_minimum_gridsearch(xtal::Crystal, molecule::Molecule{Cart}, ljff::LJForceField; + resolution::Union{Float64, Tuple{Int, Int, Int}}=1.0)::Tuple{Molecule{Frac}, Float64} # Perform an energy grid calculation on a course grid to get initial estimate. - grid = energy_grid(xtal, molecule, ljff; n_pts=n_pts) + grid = energy_grid(xtal, molecule, ljff; resolution=resolution) E_min, idx_min = findmin(grid.data) xf_minE = id_to_xf(Tuple(idx_min), grid.n_pts) # fractional coords of min diff --git a/src/eos.jl b/src/eos.jl index 4170621b8..2f53f268f 100644 --- a/src/eos.jl +++ b/src/eos.jl @@ -1,9 +1,6 @@ # Calculates the properties of a real fluid, such as the compressibility factor, fugacity, # and molar volume. -# Universal fluid constant (R). units: m³-bar/(K-mol) -const R = 8.3144598e-5 - # Data structure stating characteristics of a Peng-Robinson fluid struct PengRobinsonFluid "Peng-Robinson Fluid species. e.g. :CO2" @@ -16,14 +13,16 @@ struct PengRobinsonFluid ω::Float64 end + # Parameters in the Peng-Robinson Equation of State # T in Kelvin, P in bar -a(fluid::PengRobinsonFluid) = (0.457235 * R ^ 2 * fluid.Tc ^ 2) / fluid.Pc -b(fluid::PengRobinsonFluid) = (0.0777961 * R * fluid.Tc) / fluid.Pc +a(fluid::PengRobinsonFluid) = (0.457235 * UNIV_GAS_CONST ^ 2 * fluid.Tc ^ 2) / fluid.Pc +b(fluid::PengRobinsonFluid) = (0.0777961 * UNIV_GAS_CONST * fluid.Tc) / fluid.Pc κ(fluid::PengRobinsonFluid) = 0.37464 + (1.54226 * fluid.ω) - (0.26992 * fluid.ω ^ 2) α(κ::Float64, Tr::Float64) = (1 + κ * (1 - √Tr)) ^ 2 -A(T::Float64, P::Float64, fluid::PengRobinsonFluid) = α(κ(fluid), T / fluid.Tc) * a(fluid) * P / (R ^ 2 * T ^ 2) -B(T::Float64, P::Float64, fluid::PengRobinsonFluid) = b(fluid) * P / (R * T) +A(T::Float64, P::Float64, fluid::PengRobinsonFluid) = α(κ(fluid), T / fluid.Tc) * a(fluid) * P / (UNIV_GAS_CONST ^ 2 * T ^ 2) +B(T::Float64, P::Float64, fluid::PengRobinsonFluid) = b(fluid) * P / (UNIV_GAS_CONST * T) + # Calculates three outputs for compressibility factor using the polynomial form of # the Peng-Robinson Equation of State. Filters for only real roots and returns the @@ -44,6 +43,7 @@ function compressibility_factor(fluid::PengRobinsonFluid, T::Float64, P::Float64 return real(z_factor[id_closest_to_unity]) end + # Calculating for fugacity coefficient from an integration (bar). function calculate_ϕ(fluid::PengRobinsonFluid, T::Float64, P::Float64) z = compressibility_factor(fluid, T, P) @@ -53,11 +53,12 @@ function calculate_ϕ(fluid::PengRobinsonFluid, T::Float64, P::Float64) return exp(log_ϕ) end + """ fluid = PengRobinsonFluid(fluid) Reads in critical temperature, critical pressure, and acentric factor of the `fluid::Symbol` -from the properties .csv file `joinpath(PorousMaterials.PATH_TO_DATA, "PengRobinson_fluid_props.csv")` +from the properties .csv file `joinpath(PorousMaterials.rc[:paths][:data], "PengRobinson_fluid_props.csv")` and returns a complete `PengRobinsonFluid` data structure. **NOTE: Do not delete the last three comment lines in PengRobinson_fluid_props.csv @@ -68,10 +69,10 @@ and returns a complete `PengRobinsonFluid` data structure. - `PengRobinsonFluid::struct`: Data structure containing Peng-Robinson fluid parameters. """ function PengRobinsonFluid(fluid::Symbol) - df = CSV.read(joinpath(PATH_TO_DATA, "PengRobinson_fluid_props.csv"), DataFrame, copycols=true, comment="#") + df = CSV.read(joinpath(rc[:paths][:data], "PengRobinson_fluid_props.csv"), DataFrame, copycols=true, comment="#") filter!(row -> row[:fluid] == string(fluid), df) if nrow(df) == 0 - error(@sprintf("fluid %s properties not found in %sPengRobinson_fluid_props.csv", fluid, PATH_TO_DATA)) + error(@sprintf("fluid %s properties not found in %sPengRobinson_fluid_props.csv", fluid, rc[:paths][:data])) end Tc = df[1, Symbol("Tc(K)")] Pc = df[1, Symbol("Pc(bar)")] @@ -79,6 +80,7 @@ function PengRobinsonFluid(fluid::Symbol) return PengRobinsonFluid(fluid, Tc, Pc, ω) end + # Prints resulting values for Peng-Robinson fluid properties function Base.show(io::IO, fluid::PengRobinsonFluid) println(io, "fluid species: ", fluid.fluid) @@ -87,6 +89,7 @@ function Base.show(io::IO, fluid::PengRobinsonFluid) println(io, "\tAcenteric factor: ", fluid.ω) end + # Data structure stating characteristics of a van der Waals fluid struct VdWFluid "van der Waals Fluid species. e.g. :CO2" @@ -97,6 +100,7 @@ struct VdWFluid b::Float64 end + # Calculates the compressibility factor Z for fluids function compressibility_factor(fluid::VdWFluid, T::Float64, P::Float64) # build polynomial in ρ: D ρ³ + C ρ² + B ρ + A = 0 @@ -106,7 +110,7 @@ function compressibility_factor(fluid::VdWFluid, T::Float64, P::Float64) D = - fluid.a * fluid.b C = fluid.a - B = - (P * fluid.b + R * T) + B = - (P * fluid.b + UNIV_GAS_CONST * T) A = P # Creates polynomial in ρ the VdW cubic function @@ -119,25 +123,25 @@ function compressibility_factor(fluid::VdWFluid, T::Float64, P::Float64) # is the density corresponding to the fluid phase ρ = minimum(real_rho) # Compressibility factor - z = P / (ρ * R * T) + z = P / (ρ * UNIV_GAS_CONST * T) return z end + # Calculates for fugacity using derivation of van der Waals EOS function calculate_ϕ(fluid::VdWFluid, T::Float64, P::Float64) - z = compressibility_factor(fluid, T, P) - ρ = P / (z * R * T) - log_f = log(P) + (fluid.b - fluid.a / (R * T)) * (P / (R * T)) + log_f = log(P) + (fluid.b - fluid.a / (UNIV_GAS_CONST * T)) * (P / (UNIV_GAS_CONST * T)) # Defines the fugacity coefficient as fugacity over pressure ϕ = exp(log_f) / P return ϕ end + """ fluid = VdWFluid(fluid) Reads in van der Waals constants of the `fluid::Symbol` -from the properties .csv file `joinpath(PorousMaterials.PATH_TO_DATA, "VdW_fluid_props.csv")` +from the properties .csv file `joinpath(PorousMaterials.rc[:paths][:data], "VdW_fluid_props.csv")` and returns a complete `VdWFluid` data structure. ***NOTE: Do not delete the last three comment lines in VdW_fluid_props.csv @@ -148,16 +152,17 @@ and returns a complete `VdWFluid` data structure. - `VdWFluid::struct`: Data structure containing van der Waals constants """ function VdWFluid(fluid::Symbol) - df = CSV.read(joinpath(PATH_TO_DATA, "VdW_fluid_props.csv"), DataFrame, copycols=true, comment="#") + df = CSV.read(joinpath(rc[:paths][:data], "VdW_fluid_props.csv"), DataFrame, copycols=true, comment="#") filter!(row -> row[:fluid] == string(fluid), df) if nrow(df) == 0 - error(@sprintf("Fluid %s constants not found in %sVdW_fluidops.csv", fluid, PATH_TO_DATA)) + error(@sprintf("Fluid %s constants not found in %sVdW_fluidops.csv", fluid, rc[:paths][:data])) end a = df[1, Symbol("a(bar*m^6/mol^2)")] b = df[1, Symbol("b(m^3/mol)")] return VdWFluid(fluid, a, b) end + # Prints resulting values for van der Waals constants function Base.show(io::IO, fluid::VdWFluid) println(io, "Fluid species: ", fluid.fluid) @@ -185,7 +190,7 @@ function calculate_properties(fluid::Union{PengRobinsonFluid, VdWFluid}, T::Floa # Compressbility factor (unitless) z = compressibility_factor(fluid, T, P) # Density (mol/m^3) - ρ = P / (z * R * T) + ρ = P / (z * UNIV_GAS_CONST * T) # Molar volume (L/mol) Vm = 1000.0 / ρ # Fugacity (bar) diff --git a/src/forcefield.jl b/src/forcefield.jl index f5bb9f52e..b4f139675 100644 --- a/src/forcefield.jl +++ b/src/forcefield.jl @@ -21,8 +21,10 @@ struct LJForceField r²_cutoff::Float64 end + Base.broadcastable(ljff::LJForceField) = Ref(ljff) + """ ljforcefield = LJForceField(forcefield; r_cutoff=14.0, mixing_rules="Lorentz-Berthelot") @@ -47,7 +49,7 @@ function LJForceField(forcefield::String; r_cutoff::Float64=14.0, error(@sprintf("%s mixing rules not implemented...\n", mixing_rules)) end - forcefield_file_path = joinpath(PATH_TO_FORCEFIELDS, forcefield * ".csv") + forcefield_file_path = joinpath(rc[:paths][:forcefields], forcefield * ".csv") df = CSV.read(forcefield_file_path, DataFrame, comment="#") # from DataFrames @@ -97,6 +99,7 @@ function LJForceField(forcefield::String; r_cutoff::Float64=14.0, return ljff end + """ repfactors = replication_factors(unitcell, r_cutoff) @@ -158,6 +161,7 @@ replication_factors(unitcell::Box, ljforcefield::LJForceField) = replication_fac replication_factors(crystal::Crystal, r_cutoff::Float64) = replication_factors(crystal.box, r_cutoff) replication_factors(crystal::Crystal, ljforcefield::LJForceField) = replication_factors(crystal.box, sqrt(ljforcefield.r²_cutoff)) + """ forcefield_coverage(atoms, ljforcefield) forcefield_coverage(molecule, ljforcefield) @@ -186,6 +190,7 @@ function forcefield_coverage(atoms::Atoms, ljff::LJForceField) end forcefield_coverage(crystal::Crystal, ljff::LJForceField) = forcefield_coverage(crystal.atoms, ljff) + function Base.show(io::IO, ff::LJForceField) println(io, "Force field: ", ff.name) println(io, "Number of atoms included: ", length(ff.pure_σ)) diff --git a/src/gcmc.jl b/src/gcmc.jl index aca5a8f7e..29a2f3837 100644 --- a/src/gcmc.jl +++ b/src/gcmc.jl @@ -1,7 +1,5 @@ import Base: +, / -const KB = 1.38064852e7 # Boltmann constant (Pa-m3/K --> Pa-A3/K) - ### # Markov chain proposals ### @@ -19,12 +17,14 @@ const TRANSLATION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["translation const ROTATION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["rotation"] const REINSERTION = Dict([v => k for (k, v) in PROPOSAL_ENCODINGS])["reinsertion"] + # count proposed/accepted for each subtype mutable struct MarkovCounts n_proposed::Array{Int, 1} n_accepted::Array{Int, 1} end + ### # collecting statistics ### @@ -42,6 +42,7 @@ end GCMCstats() = GCMCstats(0, 0, 0, SystemPotentialEnergy(), SystemPotentialEnergy(), 0.0) + +(s1::GCMCstats, s2::GCMCstats) = GCMCstats(s1.n_samples + s2.n_samples, s1.n + s2.n, s1.n² + s2.n², @@ -49,6 +50,7 @@ GCMCstats() = GCMCstats(0, 0, 0, SystemPotentialEnergy(), SystemPotentialEnergy( s1.U² + s2.U², s1.Un + s2.Un) + function Base.sum(gcmc_stats::Array{GCMCstats, 1}) sum_stats = GCMCstats() for gs in gcmc_stats @@ -57,6 +59,7 @@ function Base.sum(gcmc_stats::Array{GCMCstats, 1}) return sum_stats end + function Base.print(gcmc_stats::GCMCstats) println("\t# samples: ", gcmc_stats.n_samples) println("\t⟨N⟩ (molecules) = ", gcmc_stats.n / gcmc_stats.n_samples) @@ -69,6 +72,7 @@ function Base.print(gcmc_stats::GCMCstats) println("\t⟨U⟩ (K) = ", sum(gcmc_stats.U) / gcmc_stats.n_samples) end + # Compute average and standard error of the number of molecules and potential # energy from an array of `GCMCstats`, each corresponding to statitics from an # independent block/bin during the simulation. The average from each bin is @@ -92,6 +96,7 @@ function mean_stderr_n_U(gcmc_stats::Array{GCMCstats, 1}) return avg_n, err_n, avg_U, err_U end + # TODO move this to MC helpers? but not sure if it will inline. so wait after test with @time # potential energy change after inserting/deleting/perturbing coordinates of molecules[molecule_id] @inline function potential_energy(molecule_id::Int, @@ -117,6 +122,7 @@ end return energy end + """ results, molecules = μVT_sim(xtal, molecule, temperature, pressure, ljff; molecules=Molecule[], settings=settings) @@ -331,7 +337,7 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, #### # proposal probabilities ### - mc_proposal_probabilities = [0.0 for p = 1:N_PROPOSAL_TYPES] + mc_proposal_probabilities = [0.0 for _ ∈ 1:N_PROPOSAL_TYPES] # set defaults mc_proposal_probabilities[INSERTION] = 0.35 mc_proposal_probabilities[DELETION] = mc_proposal_probabilities[INSERTION] # must be equal @@ -354,7 +360,7 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, end # initiate GCMC statistics for each block # break simulation into `N_BLOCKS` blocks to gauge convergence - gcmc_stats = [GCMCstats() for block_no = 1:N_BLOCKS] + gcmc_stats = [GCMCstats() for _ ∈ 1:N_BLOCKS] current_block = 1 # make sure the number of sample cycles is at least equal to N_BLOCKS if n_sample_cycles < N_BLOCKS @@ -373,7 +379,7 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, if show_progress_bar next!(progress_bar; showvalues=[(:cycle, outer_cycle), (:number_of_molecules, length(molecules))]) end - for inner_cycle = 1:max(20, length(molecules)) + for _ ∈ 1:max(20, length(molecules)) markov_chain_time += 1 # choose proposed move randomly; keep track of proposals @@ -387,7 +393,7 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, ljff, eparams, eikr_gh, eikr_gg) # Metropolis Hastings Acceptance for Insertion - if rand() < fugacity * xtal.box.Ω / (length(molecules) * KB * + if rand() < fugacity * xtal.box.Ω / (length(molecules) * BOLTZMANN * temperature) * exp(-sum(ΔE) / temperature) # accept the move, adjust current_energy markov_counts.n_accepted[which_move] += 1 @@ -406,7 +412,7 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, eparams, eikr_gh, eikr_gg) # Metropolis Hastings Acceptance for Deletion - if rand() < length(molecules) * KB * temperature / (fugacity * + if rand() < length(molecules) * BOLTZMANN * temperature / (fugacity * xtal.box.Ω) * exp(sum(ΔE) / temperature) # accept the deletion, delete molecule, adjust current_energy markov_counts.n_accepted[which_move] += 1 @@ -630,7 +636,6 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, results["err ⟨U_gg, electro⟩ (K)"] = err_U.gg.es results["err ⟨U⟩ (K)"] = sum(err_U) - # average N in more common units results["⟨N⟩ (molecules/unit cell)"] = avg_n / (repfactors[1] * repfactors[2] * repfactors[3]) results["err ⟨N⟩ (molecules/unit cell)"] = err_n / (repfactors[1] * repfactors[2] * repfactors[3]) @@ -659,11 +664,11 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, molecules = Cart.(molecules, xtal.box) if autosave - if ! isdir(PATH_TO_SIMS) - mkdir(PATH_TO_SIMS) + if ! isdir(rc[:paths][:sims]) + mkdir(rc[:paths][:sims]) end - save_results_filename = joinpath(PATH_TO_SIMS, + save_results_filename = joinpath(rc[:paths][:sims], μVT_output_filename(xtal, molecule, temperature, pressure, ljff, n_burn_cycles, n_sample_cycles, comment=results_filename_comment @@ -679,6 +684,7 @@ function μVT_sim(xtal::Crystal, molecule::Molecule, temperature::Float64, return results, molecules # summary of statistics and ending configuration of molecules end # μVT_sim + """ filename = μVT_output_filename(xtal, molecule, temperature, pressure, ljff, n_burn_cycles, @@ -716,6 +722,7 @@ function μVT_output_filename(xtal::Crystal, molecule::Molecule, temperature::Fl ) * comment * extension end + function print_results(results::Dict; print_title::Bool=true) if print_title # already print in GCMC tests... @@ -733,7 +740,7 @@ function print_results(results::Dict; print_title::Bool=true) println(key * ": ", results[key]) end - for (proposal_id, proposal_description) in PROPOSAL_ENCODINGS + for (_, proposal_description) in PROPOSAL_ENCODINGS total_proposals = results[@sprintf("Total # %s proposals", proposal_description)] fraction_accepted = results[@sprintf("Fraction of %s proposals accepted", proposal_description)] if total_proposals > 0 @@ -757,6 +764,7 @@ function print_results(results::Dict; print_title::Bool=true) return end + function pretty_print(xtal::Crystal, molecule::Molecule, temperature::Float64, pressure::Float64, ljff::LJForceField) print("Simulating ") printstyled("(μVT)"; color=:yellow) @@ -773,6 +781,7 @@ function pretty_print(xtal::Crystal, molecule::Molecule, temperature::Float64, p println(" force field.") end + """ results = stepwise_adsorption_isotherm(xtal, molecule, temperature, pressures, ljff; n_sample_cycles=5000, @@ -816,7 +825,7 @@ function stepwise_adsorption_isotherm(xtal::Crystal, results = Dict{String, Any}[] # push results to this array molecules = Molecule{Cart}[] # initiate with empty xtal - for (i, pressure) in enumerate(pressures) + for pressure ∈ pressures result, molecules = μVT_sim(xtal, molecule, temperature, pressure, ljff, n_burn_cycles=n_burn_cycles, @@ -837,6 +846,7 @@ function stepwise_adsorption_isotherm(xtal::Crystal, return results end + """ results = adsorption_isotherm(xtal, molecule, temperature, pressures, ljff; n_sample_cycles=5000, @@ -902,6 +912,7 @@ function adsorption_isotherm(xtal::Crystal, return results[[findall(x -> x==i, ids)[1] for i = 1:length(ids)]] end + """ dataframe = isotherm_sim_results_to_dataframe(desired_props, xtal, molecule, temperature, @@ -925,7 +936,7 @@ to locate the requested output files, this function calls [`μVT_output_filename - `n_burn_cycles::Int64`: The number of burn cycles used in this simulation - `n_sample_cycles::Int64`: The number of sample cycles used in the simulations - `where_are_jld_files::Union{String, Nothing}=nothing`: The location to the simulation files. defaults to - `PorousMaterials.PATH_TO_SIMS`. + `PorousMaterials.rc[:paths][:sims]`. - `comment::String=""`: comment appended to outputfilename # Returns @@ -963,7 +974,7 @@ function isotherm_sim_results_to_dataframe(desired_props::Array{String, 1}, where_are_jld_files::Union{String, Nothing}=nothing) # determine the location of the data files if isnothing(where_are_jld_files) - where_are_jld_files = PATH_TO_SIMS + where_are_jld_files = rc[:paths][:sims] end # prepare dataframe to populate df = DataFrame() diff --git a/src/grid.jl b/src/grid.jl index 308044cc0..b82d25397 100644 --- a/src/grid.jl +++ b/src/grid.jl @@ -17,6 +17,7 @@ struct Grid{T} origin::Array{Float64, 1} end + """ voxel_id = xf_to_id(n_pts, xf) @@ -42,6 +43,7 @@ function xf_to_id(n_pts::Tuple{Int, Int, Int}, xf::Array{Float64, 1}) return voxel_id end + """ xf = id_to_xf(voxel_id, n_pts) @@ -58,6 +60,7 @@ function id_to_xf(id::Tuple{Int, Int, Int}, n_pts::Tuple{Int, Int, Int}) return [(id[k] - 1.0) / (n_pts[k] - 1.0) for k = 1:3] end + """ update_density!(grid, molecule, species) @@ -83,6 +86,7 @@ function update_density!(grid::Grid, molecules::Array{Molecule{Frac}, 1}, specie end end + # don't include molecules outside the box. function update_density!(grid::Grid, molecules::Array{Molecule{Cart}, 1}, species::Symbol) for molecule in molecules @@ -101,6 +105,7 @@ function update_density!(grid::Grid, molecules::Array{Molecule{Cart}, 1}, specie end end + """ n_pts = required_n_pts(box, dx) @@ -117,6 +122,7 @@ function required_n_pts(box::Box, dx::Float64) return Tuple(n_pts) end + """ write_cube(grid, filename, verbose=true) @@ -126,13 +132,13 @@ The atoms of the unit cell are not printed in the .cube. Instead, use .xyz files # Arguments - `grid::Grid`: grid with associated data at each grid point. -- `filename::AbstractString`: name of .cube file to which we write the grid; this is relative to `PATH_TO_GRIDS`. +- `filename::AbstractString`: name of .cube file to which we write the grid; this is relative to `rc[:paths][:grids]`. - `verbose::Bool`: print name of file after writing. - `length_units::String`: units for length. Bohr or Angstrom. """ function write_cube(grid::Grid, filename::AbstractString; verbose::Bool=true, length_units::String="Angstrom") - if ! isdir(PATH_TO_GRIDS) - mkdir(PATH_TO_GRIDS) + if ! isdir(rc[:paths][:grids]) + mkdir(rc[:paths][:grids]) end @assert (length_units in ["Angstrom", "Bohr"]) @@ -140,7 +146,7 @@ function write_cube(grid::Grid, filename::AbstractString; verbose::Bool=true, le if ! occursin(".cube", filename) filename = filename * ".cube" end - cubefile = open(joinpath(PATH_TO_GRIDS, filename), "w") + cubefile = open(joinpath(rc[:paths][:grids], filename), "w") @printf(cubefile, "Units of data: %s\nLoop order: x, y, z\n", grid.units) @@ -171,18 +177,19 @@ function write_cube(grid::Grid, filename::AbstractString; verbose::Bool=true, le end # loop over x points close(cubefile) if verbose - println("\tSee ", joinpath(PATH_TO_GRIDS, filename)) + println("\tSee ", joinpath(rc[:paths][:grids], filename)) end return end + """ grid = read_cube(filename) Read a .cube file and return a populated `Grid` data structure. # Arguments -- `filename::AbstractString`: name of .cube file to which we write the grid; this is relative to `PATH_TO_GRIDS` +- `filename::AbstractString`: name of .cube file to which we write the grid; this is relative to `rc[:paths][:grids]` # Returns - `grid::Grid`: A grid data structure @@ -192,7 +199,7 @@ function read_cube(filename::AbstractString) filename *= ".cube" end - cubefile = open(joinpath(PATH_TO_GRIDS, filename)) + cubefile = open(joinpath(rc[:paths][:grids], filename)) # waste two lines line = readline(cubefile) @@ -204,7 +211,7 @@ function read_cube(filename::AbstractString) origin = [parse(Float64, split(line)[1 + i]) for i = 1:3] # read box info - box_lines = [readline(cubefile) for i = 1:3] + box_lines = [readline(cubefile) for _ ∈ 1:3] # number of grid pts n_pts = Tuple([parse(Int, split(box_lines[i])[1]) for i = 1:3]) @@ -242,6 +249,7 @@ function read_cube(filename::AbstractString) return Grid(box, n_pts, data, units, origin) end + """ grid = energy_grid(crystal, molecule, ljforcefield; resolution=1.0, temperature=298.0, n_rotations=750) @@ -316,8 +324,8 @@ function energy_grid(crystal::Crystal, molecule::Molecule{Cart}, ljforcefield::L ensemble_average_energy = vdw_energy(crystal, molecule, ljforcefield) else boltzmann_factor_sum = 0.0 - for r = 1:n_rotations - rotate!(molecule) + for _ ∈ 1:n_rotations + random_rotation!(molecule) energy = PotentialEnergy(0.0, 0.0) energy.vdw = vdw_energy(crystal, molecule, ljforcefield) @@ -339,12 +347,14 @@ function energy_grid(crystal::Crystal, molecule::Molecule{Cart}, ljforcefield::L return grid end + function Base.show(io::IO, grid::Grid) @printf(io, "Regular grid of %d by %d by %d points superimposed over a unit cell and associated data.\n", grid.n_pts[1], grid.n_pts[2], grid.n_pts[3]) @printf(io, "\tunits of data attribute: %s\n", grid.units) @printf(io, "\torigin: [%f, %f, %f]\n", grid.origin[1], grid.origin[2], grid.origin[3]) end + # comparing very large numbers in grid.data, so increase rtol to account # for loss of precision when writing grid.data to a cube file. function Base.isapprox(g1::Grid, g2::Grid; atol::Real=0.0, rtol::Real=atol > 0.0 ? 0.0 : sqrt(eps())) @@ -355,6 +365,7 @@ function Base.isapprox(g1::Grid, g2::Grid; atol::Real=0.0, rtol::Real=atol > 0.0 isapprox(g1.origin, g2.origin, atol=atol, rtol=rtol)) end + function _flood_fill!(grid::Grid, segmented_grid::Grid, queue_of_grid_pts::Array{Tuple{Int, Int, Int}, 1}, i::Int, j::Int, k::Int, energy_tol::Float64) @@ -392,6 +403,7 @@ function _flood_fill!(grid::Grid, segmented_grid::Grid, return nothing end + function _segment_grid(grid::Grid, energy_tol::Float64, verbose::Bool) # grid of Int's corresponding to each original grid point. # let "0" be "unsegmented" @@ -429,6 +441,7 @@ function _segment_grid(grid::Grid, energy_tol::Float64, verbose::Bool) return segmented_grid end + # this returns number of segments in a segmented grid. # it excludes the inaccessible portions -1. function _count_segments(segmented_grid::Grid) @@ -444,6 +457,7 @@ function _count_segments(segmented_grid::Grid) return nb_segments end + # struct describing directed edge describing connection between two segments of a flood- # filled grid. the direction (which unit cell face is traversed in the connection) is # also included. @@ -453,6 +467,7 @@ struct SegmentConnection direction::Tuple{Int, Int, Int} # unit cell face traversed end + function _note_connection!(segment_1::Int, segment_2::Int, connections::Array{SegmentConnection, 1}, direction::Tuple{Int, Int, Int}, verbose::Bool=true) # ignore connections between un-occupiable regions @@ -475,6 +490,7 @@ function _note_connection!(segment_1::Int, segment_2::Int, connections::Array{Se return nothing end + # obtain set of Segment Connections function _build_list_of_connections(segmented_grid::Grid; verbose::Bool=true) # initialize empty list of edges @@ -499,6 +515,7 @@ function _build_list_of_connections(segmented_grid::Grid; verbose::Bool=true) return connections end + function _translate_into_graph(segmented_grid::Grid, connections::Array{SegmentConnection, 1}) nb_segments = _count_segments(segmented_grid) @@ -534,6 +551,7 @@ function _translate_into_graph(segmented_grid::Grid, connections::Array{SegmentC return graph, vertex_to_direction end + function _classify_segments(segmented_grid::Grid, graph::SimpleDiGraph{Int64}, vertex_to_direction::Dict{Int, Union{Nothing, Tuple{Int, Int, Int}}}, verbose::Bool=true) @@ -550,8 +568,8 @@ function _classify_segments(segmented_grid::Grid, graph::SimpleDiGraph{Int64}, # all edges involved in a cycle are connected together and are accessible channels for cyc in all_cycles - @assert vertex_to_direction[cyc[1]] == nothing "shld start with segment." - @assert vertex_to_direction[cyc[end]] != nothing "shld end with direction vertex." + @assert isnothing(vertex_to_direction[cyc[1]]) "should start with segment." + @assert !isnothing(vertex_to_direction[cyc[end]]) "should end with direction vertex." # travel through cycle, keep track of which unit cell we're in # start in home unit cell unit_cell = (0, 0, 0) @@ -559,7 +577,7 @@ function _classify_segments(segmented_grid::Grid, graph::SimpleDiGraph{Int64}, for s in cyc # if an intermediate unit cell boundary vertex, then keep track of unit cell # we're in by looking at boundary we traversed. - if vertex_to_direction[s] != nothing + if !isnothing(vertex_to_direction[s]) unit_cell = unit_cell .+ vertex_to_direction[s] end end @@ -571,7 +589,7 @@ function _classify_segments(segmented_grid::Grid, graph::SimpleDiGraph{Int64}, println("\t...found a cycle of accessible segments!") end for s in cyc - if vertex_to_direction[s] == nothing # i.e. if this is not a direciton vertex + if isnothing(vertex_to_direction[s]) # i.e. if this is not a direciton vertex segment_classifiction[s] = 1 end end @@ -596,6 +614,7 @@ function _classify_segments(segmented_grid::Grid, graph::SimpleDiGraph{Int64}, return segment_classifiction end + function _assign_inaccessible_pockets_minus_one!(segmented_grid::Grid, segment_classifiction::Array{Int, 1}; verbose::Bool=true) for s = 1:length(segment_classifiction) @@ -612,6 +631,7 @@ function _assign_inaccessible_pockets_minus_one!(segmented_grid::Grid, end end + """ accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal, probe_molecule, ljforcefield; resolution::Union{Float64, Tuple{Int, Int, Int}}=1.0, energy_tol=10.0, energy_unit=:kJ_mol, @@ -739,11 +759,13 @@ function compute_accessibility_grid(crystal::Crystal, probe::Molecule, forcefiel return accessibility_grid, nb_segments_blocked, porosity end + # return ID that is the nearest neighbor. function _arg_nearest_neighbor(n_pts::Tuple{Int, Int, Int}, xf::Array{Float64, 1}) return 1 .+ round.(Int, xf .* (n_pts .- 1)) end + function _apply_pbc_to_index!(id::Array{Int, 1}, n_pts::Tuple{Int, Int, Int}) for k = 1:3 if id[k] == 0 @@ -756,6 +778,7 @@ function _apply_pbc_to_index!(id::Array{Int, 1}, n_pts::Tuple{Int, Int, Int}) return nothing end + """ accessible(accessibility_grid, xf) accessible(accessibility_grid, xf, repfactors) @@ -797,3 +820,5 @@ end function accessible(accessibility_grid::Grid{Bool}, xf::Array{Float64, 1}, repfactors::Tuple{Int, Int, Int}) return accessible(accessibility_grid, mod.(xf .* repfactors, 1.0)) end + +origin(T::DataType) = T([0.0, 0.0, 0.0]) \ No newline at end of file diff --git a/src/henry.jl b/src/henry.jl index c81a208d4..dced5d195 100644 --- a/src/henry.jl +++ b/src/henry.jl @@ -1,8 +1,4 @@ -const UNIVERSAL_GAS_CONSTANT = 8.3144598e-5 # m³-bar/(K-mol) -const K_to_KJ_PER_MOL = 8.3144598 / 1000.0 # kJ/(mol-K) - """ - ``` result = henry_coefficient(crystal, molecule, temperature, ljforcefield, nb_insertions=1e6, verbose=true, ewald_precision=1e-6, @@ -29,7 +25,7 @@ average, per unit cell volume (ų) - `verbose::Bool`: whether or not to print off information during the simulation. - `ewald_precision::Float64`: desired precision for Ewald summations; used to determine the replication factors in reciprocal space. -- `autosave::Bool`: save results file as a .jld in PATH_TO_DATA * `sims` +- `autosave::Bool`: save results file as a .jld in rc[:paths][:data] * `sims` - `filename_comment::AbstractString`: An optional comment that will be appended to the name of the saved file. # Returns @@ -58,7 +54,7 @@ function henry_coefficient(crystal::Crystal, molecule::Molecule, temperature::Fl printstyled(insertions_per_volume; color=:green) println(" insertions per ų.") - if accessibility_grid != nothing + if !isnothing(accessibility_grid) println("Using provided accessibility grid to block inaccessible pockets") if ! isapprox(accessibility_grid.box, crystal.box) error(@sprintf("accessibility grid box does not match box of %s.\n", @@ -121,7 +117,7 @@ function henry_coefficient(crystal::Crystal, molecule::Molecule, temperature::Fl # Kₕ = β Σ e ^(βUᵢ) / nb_insertions_per_block # (these N_BLOCKS-long arrays) average_energies = wtd_energy_sums ./ boltzmann_factor_sums # K - henry_coefficients = boltzmann_factor_sums / (UNIVERSAL_GAS_CONSTANT * temperature * nb_insertions_per_block) # mol/(m³-bar) + henry_coefficients = boltzmann_factor_sums / (UNIV_GAS_CONST * temperature * nb_insertions_per_block) # mol/(m³-bar) if verbose for b = 1:N_BLOCKS @@ -152,13 +148,13 @@ function henry_coefficient(crystal::Crystal, molecule::Molecule, temperature::Fl results["⟨U, es⟩ (K)"] = mean([average_energies[b].es for b = 1:N_BLOCKS]) results["⟨U⟩ (K)"] = results["⟨U, vdw⟩ (K)"] + results["⟨U, es⟩ (K)"] - results["⟨U⟩ (kJ/mol)"] = results["⟨U⟩ (K)"] * K_to_KJ_PER_MOL - results["⟨U, vdw⟩ (kJ/mol)"] = results["⟨U, vdw⟩ (K)"] * K_to_KJ_PER_MOL - results["err ⟨U, vdw⟩ (kJ/mol)"] = err_energy.vdw * K_to_KJ_PER_MOL - results["⟨U, es⟩ (kJ/mol)"] = results["⟨U, es⟩ (K)"] * K_to_KJ_PER_MOL - results["err ⟨U, es⟩ (kJ/mol)"] = err_energy.es * K_to_KJ_PER_MOL - results["Qst (kJ/mol)"] = -results["⟨U⟩ (kJ/mol)"] + temperature * K_to_KJ_PER_MOL - results["err Qst (kJ/mol)"] = sum(err_energy) * K_to_KJ_PER_MOL + results["⟨U⟩ (kJ/mol)"] = results["⟨U⟩ (K)"] * K_TO_KJ_PER_MOL + results["⟨U, vdw⟩ (kJ/mol)"] = results["⟨U, vdw⟩ (K)"] * K_TO_KJ_PER_MOL + results["err ⟨U, vdw⟩ (kJ/mol)"] = err_energy.vdw * K_TO_KJ_PER_MOL + results["⟨U, es⟩ (kJ/mol)"] = results["⟨U, es⟩ (K)"] * K_TO_KJ_PER_MOL + results["err ⟨U, es⟩ (kJ/mol)"] = err_energy.es * K_TO_KJ_PER_MOL + results["Qst (kJ/mol)"] = -results["⟨U⟩ (kJ/mol)"] + temperature * K_TO_KJ_PER_MOL + results["err Qst (kJ/mol)"] = sum(err_energy) * K_TO_KJ_PER_MOL results["elapsed time (min)"] = elapsed_time / 60 diff --git a/src/matter.jl b/src/matter.jl deleted file mode 100644 index 7f80f6afd..000000000 --- a/src/matter.jl +++ /dev/null @@ -1,217 +0,0 @@ -### -# coordinates: Fractional and Cartesian -### -""" -abstract type for coordinates. -""" -abstract type Coords end - -#Base.IndexStyle(::Type{<:Coords}) = IndexLinear() - -""" -fractional coordinates, a subtype of `Coords`. - -construct by passing an `Array{Float64, 2}` whose columns are the coordinates. - -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 - xf::Array{Float64, 2} -end -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)) - -""" -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 - x::Array{Float64, 2} -end -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 -Frac(xf::Array{Float64, 1}) = Frac(reshape(xf, (3, 1))) -Cart(x::Array{Float64, 1}) = Cart(reshape(x, (3, 1))) - -Base.getindex(coords::Frac, ids) = Frac(coords.xf[:, ids]) -Base.getindex(coords::Cart, ids) = Cart(coords.x[:, ids]) - -Base.length(coords::Cart) = size(coords.x, 2) -Base.length(coords::Frac) = size(coords.xf, 2) -Base.lastindex(coords::Coords) = length(coords) - -Base.setindex!(coords::Frac, val::Array{Float64, 1}, ids) = (coords.xf[:, ids] = val) -Base.setindex!(coords::Cart, val::Array{Float64, 1}, ids) = (coords.x[:, ids] = val) - -Base.size(coords::Cart) = size(coords.x) -Base.size(coords::Frac) = size(coords.xf) - -origin(T::DataType) = T([0.0, 0.0, 0.0]) - -""" - translate_by!(coords, dx) - translate_by!(coords, dx, box) - translate_by!(molecule, dx) - translate_by!(molecule, dx, box) - -translate `coords` by the vector `dx`. that is, add the vector `dx`. - -this works for any combination of `Frac` and `Cart` coords. - -modifies coordinates in place. - -`box` is needed when mixing `Frac` and `Cart` coords. - -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 - -""" - wrap!(f::Frac) - wrap!(crystal::Crystal) - -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 - -### -# Atoms -### -# T is either Cart or Frac -""" -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` - n::Int # how many atoms? - species::Array{Symbol, 1} # list of species - coords::T # coordinates -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` - n::Int # how many atoms? - species::Array{Symbol, 1} # list of species - coords::T # coordinates -end - -function Atoms(species::Array{Symbol, 1}, coords::Coords) - @assert length(species) == size(coords)[2] - @assert size(coords)[1] == 3 - return Atoms(length(species), species, 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 - -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)) - -Base.getindex(atoms::Atoms, ids) = Atoms(atoms.species[ids], atoms.coords[ids]) -Base.lastindex(atoms::Atoms) = atoms.n - -### -# point charges -### -""" -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` - n::Int - q::Array{Float64, 1} - coords::T -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` - n::Int - q::Array{Float64, 1} - coords::T -end -function Charges(q::Array{Float64, 1}, coords::Coords) - @assert length(q) == size(coords)[2] - @assert size(coords)[1] == 3 - return Charges(length(q), q, 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 - -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)) - -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) - nc = net_charge(molecule) - -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) - if charges.n == 0 - return 0.0 - else - return sum(charges.q) - end -end - -""" - neutral(charges, tol) # true or false. default tol = 1e-5 - neutral(crystal, tol) # true or false. default tol = 1e-5 - -determine if a set of `charges::Charges` (`charges.q`) sum to an absolute value less than `tol::Float64`. if `crystal::Crystal` is passed, the function looks at the `crystal.charges`. i.e. determine the absolute value of the net charge is less than `tol`. -""" -neutral(charges::Charges, tol::Float64=1e-5) = abs(net_charge(charges)) < tol diff --git a/src/misc.jl b/src/misc.jl deleted file mode 100644 index 397bf7ee7..000000000 --- a/src/misc.jl +++ /dev/null @@ -1,108 +0,0 @@ -function add_extension(filename::String, extension::String) - if ! occursin(extension, filename) - filename *= extension - end - return filename -end - -""" - atoms = read_xyz(filename) - -read a list of atomic species and their corresponding coordinates from an .xyz file. - -# Arguments -- `filename::AbstractString`: The filename of the .xyz file - -# Returns -- `atoms::Atoms{Cart}`: the set of atoms read from the .xyz file. -""" -function read_xyz(filename::AbstractString) - f = open(filename) - lines = readlines(f) - if length(lines) == 0 - return Symbol[], Float64[] - end - n = parse(Int, lines[1]) # get number of atoms - species = Symbol[] - x = zeros(Float64, 3, n) - for i = 1:n - push!(species, Symbol(split(lines[i + 2])[1])) - for j = 1:3 - x[j, i] = parse(Float64, split(lines[i + 2])[1 + j]) - end - end - close(f) - return Atoms(species, Cart(x)) -end - -""" - write_xyz(atoms, filename; comment="") - write_xyz(crystal; comment="", center_at_origin=false) - write_xyz(molecules, box, filename; comment="") # fractional - write_xyz(molecules, box, filename; comment="") # Cartesian - -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. -""" -function write_xyz(atoms::Atoms{Cart}, filename::AbstractString; comment::AbstractString="") - filename = add_extension(filename, ".xyz") - - xyzfile = open(filename, "w") - @printf(xyzfile, "%d\n%s\n", atoms.n, comment) - for i = 1:atoms.n - @printf(xyzfile, "%s\t%.4f\t%.4f\t%.4f\n", atoms.species[i], - atoms.coords.x[1, i], atoms.coords.x[2, i], atoms.coords.x[3, i]) - end - close(xyzfile) - return nothing -end - -""" - atom_colors = read_cpk_colors() - -Read in CPK color scheme for atoms. Return `atom_colors::Dict{Symbol, Tuple{Int, Int, Int}}` such that -`atom_colors[":C"]` gives RGB code for carbon as a tuple, `(144, 144, 144)`. -https://en.wikipedia.org/wiki/CPK_coloring - -# Returns -- `atom_colors::Dict{Symbol, Tuple{Int, Int, Int}}`: A dictionary linking an element symbol to its' corresponding CPK color in RGB -""" -function read_cpk_colors() - atom_colors = Dict{Symbol, Tuple{Int, Int, Int}}() - df_colors = CSV.read(joinpath(PATH_TO_DATA, "cpk_atom_colors.csv"), DataFrame) - for row in eachrow(df_colors) - atom_colors[Symbol(row[:atom])] = (row[:R], row[:G], row[:B]) - end - return atom_colors -end - -""" - atomic_masses = read_atomic_masses() - -Read the `data/atomicmasses.csv` file to construct a dictionary of atoms and their atomic -masses in amu. - -# Returns -- `atomic_masses::Dict{Symbol, Float64}`: A dictionary containing the atomic masses of each atom stored in `data/atomicmasses.csv` -""" -function read_atomic_masses() - if ! isfile(joinpath(PATH_TO_DATA, "atomicmasses.csv")) - error("Cannot find atomicmasses.csv file in your data folder\n") - end - - df_am = CSV.read(joinpath(PATH_TO_DATA, "atomicmasses.csv"), DataFrame) - - atomic_masses = Dict{Symbol, Float64}() - - for row in eachrow(df_am) - atomic_masses[Symbol(row[:atom])] = row[:mass] - end - - return atomic_masses -end diff --git a/src/molecule.jl b/src/molecule.jl index caa540190..67386d8dc 100644 --- a/src/molecule.jl +++ b/src/molecule.jl @@ -1,7 +1,3 @@ -function _define_atomic_mass() - global ATOMIC_MASS = read_atomic_masses() # for center-of-mass calcs -end - """ Data structure for a molecule/adsorbate. @@ -18,31 +14,30 @@ struct Molecule{T} # T = Frac or Cart com::T # center of mass end + function Base.isapprox(m1::Molecule, m2::Molecule) return (m1.species == m2.species) && isapprox(m1.com, m2.com) && isapprox(m1.atoms, m2.atoms) && isapprox(m1.charges, m2.charges) end + function center_of_mass(molecule::Molecule{Cart}) - if ! @isdefined ATOMIC_MASS - _define_atomic_mass() - end total_mass = 0.0 # total mass x_com = [0.0, 0.0, 0.0] # center of mass for a = 1:molecule.atoms.n - total_mass += ATOMIC_MASS[molecule.atoms.species[a]] - x_com += ATOMIC_MASS[molecule.atoms.species[a]] * molecule.atoms.coords.x[:, a] + total_mass += rc[:atomic_masses][molecule.atoms.species[a]] + x_com += rc[:atomic_masses][molecule.atoms.species[a]] * molecule.atoms.coords.x[:, a] end return Cart(x_com / total_mass) end - + + """ molecule = Molecule(species, check_neutrality=true) -construt a `Molecule` from input files read from: -`joinpath(PATH_TO_MOLECULES, species)`. +construt a `Molecule` from input files read from `joinpath(rc[:paths][:molecules], species)` -center of mass assigned using atomic masses from `read_atomic_masses()`. +center of mass assigned using atomic masses from `rc[:atomic_masses]`. # Arguments - `species::String`: name of the molecule @@ -55,7 +50,7 @@ function Molecule(species::String; check_neutrality::Bool=true) ### # Read in Lennard Jones spheres ### - df_lj = CSV.read(joinpath(PATH_TO_MOLECULES, species, "atoms.csv"), DataFrame) + df_lj = CSV.read(joinpath(rc[:paths][:molecules], species, "atoms.csv"), DataFrame) atoms = Atoms{Cart}(nrow(df_lj)) # pre-allocate atoms for (a, row) in enumerate(eachrow(df_lj)) @@ -66,7 +61,7 @@ function Molecule(species::String; check_neutrality::Bool=true) ### # Read in point charges ### - df_c = CSV.read(joinpath(PATH_TO_MOLECULES, species, "charges.csv"), DataFrame) + df_c = CSV.read(joinpath(rc[:paths][:molecules], species, "charges.csv"), DataFrame) charges = Charges{Cart}(nrow(df_c)) # pre-allocate charges for (c, row) in enumerate(eachrow(df_c)) @@ -88,9 +83,12 @@ function Molecule(species::String; check_neutrality::Bool=true) return molecule end + # documented in matter.jl +import Xtals.net_charge net_charge(molecule::Molecule) = net_charge(molecule.charges) + # convert between fractional and cartesian coords Frac(molecule::Molecule{Cart}, box::Box) = Molecule(molecule.species, Frac(molecule.atoms, box), @@ -103,6 +101,7 @@ Cart(molecule::Molecule{Frac}, box::Box) = Molecule(molecule.species, Cart(molecule.com, box) ) + """ write_xyz(box, molecules, xyz_file) @@ -110,12 +109,10 @@ Writes the coordinates of all atoms in molecules to the given xyz_file file obje passing a file object around is faster for simulation because it can be opened once at the beginning of the simulation and closed at the end. -This writes the coordinates of the molecules in cartesian coordinates, so the -box is needed for the conversion. +This writes the coordinates of the molecules in cartesian coordinates, so the box is needed for the conversion. # Arguments -- `box::Box`: The box the molecules are in, to convert molecule positions - to cartesian coordinates +- `box::Box`: The box the molecules are in, to convert molecule positions to cartesian coordinates - `molecules::Array{Molecule{Frac}, 1}`: The array of molecules to be written to the file - `xyz_file::IOStream`: The open 'write' file stream the data will be saved to """ @@ -130,7 +127,9 @@ function write_xyz(box::Box, molecules::Array{Molecule{Frac}, 1}, xyz_file::IOSt end end + # documented in matter.jl +import Xtals.translate_by! function translate_by!(molecule::Molecule{Cart}, dx::Cart) translate_by!(molecule.atoms.coords, dx) translate_by!(molecule.charges.coords, dx) @@ -151,6 +150,7 @@ function translate_by!(molecule::Molecule{Frac}, dx::Cart, box::Box) translate_by!(molecule, Frac(dx, box)) end + """ translate_to!(molecule, xf) translate_to!(molecule, x) @@ -174,6 +174,7 @@ translate_to!(molecule::Molecule{Frac}, x::Cart, box::Box) = translate_to!(molec translate_to!(molecule::Molecule{Cart}, xf::Frac, box::Box) = translate_to!(molecule, Cart(xf, box)) + function Base.show(io::IO, molecule::Molecule) println(io, "Molecule species: ", molecule.species) println(io, "Center of mass (fractional coords): ", molecule.com) @@ -207,6 +208,7 @@ function Base.show(io::IO, molecule::Molecule) end end + """ r = random_rotation_matrix() # rotation matrix in cartesian coords @@ -235,6 +237,7 @@ end random_rotation_matrix(box::Box) = box.c_to_f * random_rotation_matrix() * box.f_to_c + """ random_rotation!(molecule{Frac}, box) random_rotation!(molecule{Cart}) @@ -269,10 +272,13 @@ function random_rotation!(molecule::Molecule{Frac}, box::Box) return nothing end + # based on center of mass +import Xtals.inside inside(molecule::Molecule{Cart}, box::Box) = inside(molecule.com, box) inside(molecule::Molecule{Frac}) = inside(molecule.com) + # docstring in Misc.jl function write_xyz(molecules::Array{Molecule{Cart}, 1}, filename::AbstractString; comment::AbstractString="") @@ -297,12 +303,16 @@ function write_xyz(molecules::Array{Molecule{Frac}, 1}, box::Box, filename::Abst write_xyz(molecules, filename, comment=comment) # above end + # documented in crystal.jl +import Xtals.has_charges has_charges(molecule::Molecule) = molecule.charges.n > 0 + # documented in forcefield.jl forcefield_coverage(molecule::Molecule, ljff::LJForceField) = forcefield_coverage(molecule.atoms, ljff) + """ molecule = ion(q, coords) Facilitate constructing a point charge by constructing a molecule: @@ -319,6 +329,7 @@ function ion(q::Float64, coords::Frac) return Molecule(:ion, Atoms{Frac}(0), Charges(q, coords), coords) end + """ is_distorted = distortion(molecule, ref_molecule, box; atol=1e-12, throw_warning=true) @@ -345,21 +356,5 @@ function distortion(molecule::Molecule{Frac}, ref_molecule::Molecule{Frac}, box: pairwise_distances(ref_molecule.charges.coords, box, false), atol=atol) return true end - # loop over all pairs - # for i = 1:molecule.atoms.n - # for j = (i+1):molecule.atoms.n - # # molecule - # dxf = molecule.atoms.coords.xf[:, i] - molecule.atoms.coords.xf[:, j] - # dx = box.f_to_c * dxf - # r = norm(dx) - # # ref molecule - # dxf_ref = ref_molecule.atoms.coords.xf[:, i] - ref_molecule.atoms.coords.xf[:, j] - # dx_ref = box.f_to_c * dxf_ref - # r_ref = norm(dx_ref) - # if ! isapprox(r, r_ref, atol=atol) - # return true - # end - # end - # end return false # if made it this far... end diff --git a/test/assert_p1_symmetry.jl b/test/assert_p1_symmetry.jl deleted file mode 100644 index 86723b5f0..000000000 --- a/test/assert_p1_symmetry.jl +++ /dev/null @@ -1,44 +0,0 @@ -module Simulation_Rules_Test - -using PorousMaterials -using Test - -# Test set for making sure simulations meet certain rules -# i.e. frameworks musts be in P1 symmetry to be used in GCMC or Henry -# coefficient test -@testset "Simulation Rules Tests" begin - non_P1_framework = Crystal("symmetry_test_structure.cif", convert_to_p1=false) - molecule = Molecule("CO2") - ljff = LJForceField("UFF") - temp = 298.0 - pressure = 5.0 - - # make sure that the framework is not in P1 before running tests - - # assertion errors should be thrown when trying to run gcmc simulations - # with non-P1 frameworks - @test_throws ErrorException μVT_sim(non_P1_framework, molecule, - temp, pressure, ljff) - pressures = [1.0, 5.0, 10.0, 15.0, 20.0] - @test_throws ErrorException adsorption_isotherm(non_P1_framework, molecule, - temp, pressures, ljff) - @test_throws ErrorException stepwise_adsorption_isotherm(non_P1_framework, - molecule, temp, pressures, ljff) - - # assertion error should be thrown when attempting to run a henry - # coefficient with a non-P1 framework - @test_throws ErrorException henry_coefficient(non_P1_framework, molecule, - temp, ljff) - - # Test that an assertion is thrown when a non-P1 structure is passed into - # energy_grid(), and compute_accessibility_grid() - @test_throws ErrorException energy_grid(non_P1_framework, molecule, ljff) - @test_throws ErrorException compute_accessibility_grid(non_P1_framework, molecule, ljff) - - # Test that an assertion is thrown when trying to replicate a non-P1 - # structure - @test_throws ErrorException replicate(non_P1_framework, (2, 2, 2)) - @test_throws ErrorException replicate(non_P1_framework, (1, 1, 1)) - -end -end diff --git a/test/bonds.jl b/test/bonds.jl deleted file mode 100644 index 9106aea76..000000000 --- a/test/bonds.jl +++ /dev/null @@ -1,53 +0,0 @@ -module Bonds_test - -using PorousMaterials -using LightGraphs -using Test - -function visual_check(xtal::String) - c = Crystal(xtal) - c = replicate(c, (2, 2, 2)) - strip_numbers_from_atom_labels!(c) - infer_geometry_based_bonds!(c, true) # must - write_xyz(c) - write_bond_information(c) - @info c.name * " see .vtk and .xyz to visually check bonds" -end - -@testset "Bond 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, :*, 0.4, 1.2), - BondingRule(:Ni, :O, 0.4, 2.5), - BondingRule(:Ni, :N, 0.4, 2.5), - BondingRule(:*, :*, 0.4, 1.9)] - infer_bonds!(c, true, bonding_rules) - g1 = deepcopy(c.bonds) - conn_comps = connected_components(c.bonds) - @test length(conn_comps) == 2 # interpenetrated - c_red = getindex(c, conn_comps[1]) - c_blue = getindex(c, conn_comps[2]) - @test ne(c_red.bonds) + ne(c_blue.bonds) == ne(c.bonds) - remove_bonds!(c) - infer_geometry_based_bonds!(c, true) - @test length(conn_comps) == 2 # interpenetrated - @test g1 == c.bonds # consistency between two different bonding schemes - - # FIQCEN bonding - c = Crystal("FIQCEN_clean.cif") - strip_numbers_from_atom_labels!(c) - infer_geometry_based_bonds!(c, true) - @test length(connected_components(c.bonds)) == 1 # not interpenetrated - @test c.atoms.species[neighbors(c.bonds, 1)] == [:Cu, :O, :O, :O, :O] - visual_check("FIQCEN_clean.cif") - # reduce covalant radii to see Cu-O bond disappear - cov_radii = cordero_covalent_atomic_radii() - cov_radii[:Cu] = 1.125 - remove_bonds!(c) - @test ne(c.bonds) == 0 - infer_geometry_based_bonds!(c, true, covalent_radii=cov_radii) - @test c.atoms.species[neighbors(c.bonds, 1)] == [:O, :O, :O, :O] -end -end diff --git a/test/box.jl b/test/box.jl deleted file mode 100644 index 58f003ae8..000000000 --- a/test/box.jl +++ /dev/null @@ -1,97 +0,0 @@ -module Box_Test - -using PorousMaterials -using LinearAlgebra -using Test - -@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 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) - @test box.Ω ≈ det(box.f_to_c) - - # converters - box = Box(11.6, 5.5, 22.9, 90.0 * π / 180, 100.8 * π / 180.0, 92.0 * π / 180.0) - f = Frac([0.05, 0.9, 0.99]) - c = Cart(f, box) - f_conv = Frac(c, box) - @test isapprox(c.x, box.f_to_c * f.xf) - @test isapprox(f, f_conv) - - box = unit_cube() - @test isapprox(Box(1.0, 1.0, 1.0), unit_cube()) # alt. constructor, right angles assumed - @test box.Ω ≈ 1.0 - @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] - ) - atoms = Atoms([:C, :O], f) - atoms_c = Cart(atoms, box) - @test atoms_c.species == atoms.species - @test isapprox(atoms_c.coords.x, box.f_to_c * atoms.coords.xf) - @test isapprox(atoms, Frac(atoms_c, box)) - charges = Charges([1.0, -1.0], f) - charges_c = Cart(charges, box) - @test isapprox(charges_c.q, charges.q) - @test isapprox(charges_c.coords.x, box.f_to_c * charges.coords.xf) - @test isapprox(charges, Frac(charges_c, box)) - - # 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] - ) - @test inside(f) - f.xf[2, 2] = -0.1 - @test ! inside(f) - f.xf[2, 2] = 1.2 - @test ! inside(f) - - c = Cart([0.1 0.4; - 4.5 13.0; - 22. 10.1] - ) - @test inside(c, box) - c.x[2, 2] = -0.1 - @test !inside(c, box) - c.x[2, 2] = 20.1 - @test !inside(c, box) - - box = Box(11.6, 5.5, 22.9, 90.0 * π / 180, 100.8 * π / 180.0, 90.0 * π / 180.0) - f = Frac([0.1, 0.8, 0.9]) - @test inside(f) - @test inside(Cart(f, box), box) - - # translate_by! - box = Box(1.0, 10.0, 40.0) - 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] - ) - ) - - 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]) - ) -end -end diff --git a/test/crystal.jl b/test/crystal.jl deleted file mode 100644 index db8c8d263..000000000 --- a/test/crystal.jl +++ /dev/null @@ -1,317 +0,0 @@ -module Crystal_Test - -using PorousMaterials -using LinearAlgebra -using Test - -# 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) - ac = ac1 + ac2 - ac_dm = remove_duplicates(ac, box, true) - return isapprox(ac_dm, ac1) -end - -@testset "Crystal Tests" begin - # cif reader - xtal = Crystal("test_structure2.cif") - 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 - - # assign charges function - xtal = Crystal("test_structure3.cif") - @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 chemical_formula(xtal) == Dict(:Ca => 1, :O => 1) - @test molecular_weight(xtal) ≈ 15.9994 + 40.078 - - # overlap checker, charge neutrality checker, wrap coords checker - @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 chemical_formula(xtal) == Dict(:Ca => 2, :O => 1, :C => 1) - - # xyz writer - xtal = Crystal("SBMOF-1.cif") - 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 - - # test .cif writer; write, read in, assert equal - crystal = Crystal("SBMOF-1.cif") - rewrite_filename = "rewritten.cif" - write_cif(crystal, joinpath(PorousMaterials.PATH_TO_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(PorousMaterials.PATH_TO_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) - # rm(rewrite_filename) # clean up. - rewrite_filename = "rewritten.cif" - write_cif(crystal, joinpath(PorousMaterials.PATH_TO_CRYSTALS, rewrite_filename), number_atoms=false) - crystal_reloaded = Crystal(rewrite_filename) - @test isapprox(crystal, crystal_reloaded, atol=0.0001) - - ### 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 = Crystal("symmetry_test_structure.cif", wrap_coords=true) # default - @test all(xtal.atoms.coords.xf .< 1.0) && all(xtal.atoms.coords.xf .> 0.0) - - non_P1_crystal = Crystal("symmetry_test_structure.cif") # 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 = Crystal("symmetry_test_structure_cartn.cif") # 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) - - # 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") - # 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 - - 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) - - # 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) # TODO remove - 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 = Crystal("symmetry_test_structure.cif", convert_to_p1=true) - strip_numbers_from_atom_labels!(non_P1_xtal) - non_P1_xtal_cartesian = Crystal("symmetry_test_structure_cartn.cif", convert_to_p1=true) - 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) - - # test .cssr reader too; test_structure2.{cif,cssr} designed to be the same. - xtal_cssr = Crystal("test_structure2.cssr") - strip_numbers_from_atom_labels!(xtal_cssr) - xtal_cif = Crystal("test_structure2.cif") - strip_numbers_from_atom_labels!(xtal_cif) - @test isapprox(xtal_cif, xtal_cssr) - - # sanity checks on replicate crystal - sbmof = Crystal("SBMOF-1.cif") - 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 = replication_factors(sbmof.box, 14.0) - replicated_sbmof = replicate(sbmof, repfactors) - @test 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 chemical_formula(sbmof) == chemical_formula(replicated_sbmof) - @test isapprox(crystal_density(sbmof), crystal_density(replicated_sbmof), atol=1e-7) - - sbmof1 = Crystal("SBMOF-1.cif") - rbox = replicate(sbmof1.box, (2, 3, 4)) - @test rbox.Ω ≈ sbmof1.box.Ω * 2 * 3 * 4 - @test all(rbox.c_to_f * sbmof1.box.f_to_c * [1.0, 1.0, 1.0] .≈ [1/2, 1/3, 1/4]) - - # - # - # # test symmetry_rules - # # define default symmetry_rules - # symmetry_rules = [Array{AbstractString, 2}(undef, 3, 0) ["x", "y", "z"]] - # other_symmetry_rules = [Array{AbstractString, 2}(undef, 3, 0) ["y + z", "x + z", "x + y"]] - # symmetry_rules_two = [Array{AbstractString, 2}(undef, 3, 0) ["x" "y + z"; - # "y" "x + z"; - # "z" "x + y"]] - # symmetry_rules_two_cpy = deepcopy(symmetry_rules_two) - # @test ! is_symmetry_equal(symmetry_rules, symmetry_rules_two) - # @test ! is_symmetry_equal(symmetry_rules, other_symmetry_rules) - # @test is_symmetry_equal(symmetry_rules, symmetry_rules) - # @test is_symmetry_equal(symmetry_rules_two, symmetry_rules_two_cpy) - # - # 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 + sbmof # 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 - # TODO test bonds too. - # addition_bonding_rules = [BondingRule(:a, :b, 4.5, 5.3), - # BondingRule(:c, :d, 4.5, 5.3)] - # @test is_bonded(f1, 1, 2, [BondingRule(:a, :b, 1.0, 5.5)]; include_bonds_across_periodic_boundaries=false) - # @test ! is_bonded(f2, 1, 2, [BondingRule(:c, :d, 1.0, 4.5)]; include_bonds_across_periodic_boundaries=false) - # infer_bonds!(f1, false, addition_bonding_rules) - # infer_bonds!(f2, false, addition_bonding_rules) - # @test ! compare_bonds_in_crystal(f1 + f2, f3) - # infer_bonds!(f3, false, addition_bonding_rules) - # @test compare_bonds_in_crystal(f1 + f2, f3) - - # 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) - - # more xtal tests - sbmof1 = Crystal("SBMOF-1.cif") - @test ! has_charges(sbmof1) - @test isapprox(sbmof1.box.reciprocal_lattice, 2 * π * inv(sbmof1.box.f_to_c)) - @test sbmof1.box.Ω ≈ det(sbmof1.box.f_to_c) # sneak in crystal test - @test isapprox(crystal_density(sbmof1), 1570.4, atol=0.5) # kg/m3 - - # indexing - sbmof1_sub = sbmof1[10:15] - @test sbmof1_sub.atoms.n == 6 - @test sbmof1_sub.charges.n == 0 - @test sbmof1_sub.atoms.species == sbmof1.atoms.species[10:15] - @test sbmof1_sub.atoms.coords.xf == sbmof1.atoms.coords.xf[:, 10:15] - - # including zero charges or not, when reading in a .cif `include_zero_charges` flag to Crystal constructor - frame1 = Crystal("ATIBOU01_clean.cif") # 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"; include_zero_charges=true) - @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, PorousMaterials.NET_CHARGE_TOL) # is neutral... - write_cif(crystal, joinpath(PorousMaterials.PATH_TO_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(PorousMaterials.PATH_TO_CRYSTALS, "high_precision_charges_test_not_neutral.cif")) - @test ! neutral(Crystal("high_precision_charges_test_not_neutral.cif", check_neutrality=false), PorousMaterials.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 - # 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") - @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 = Crystal("SBMOF-1.cif") - xtal2 = Crystal("SBMOF-1.cif") - @test xtal1.atoms.species == xtal2.atoms.species - - @test_throws KeyError Crystal("SBMOF-1.cif", species_col=["bogus_column"]) - - # overlap tests - @test_throws ErrorException Crystal("SBMOF-1_overlap.cif") - xtal = Crystal("SBMOF-1_overlap.cif", remove_duplicates=true) - @test isapprox(xtal, Crystal("SBMOF-1.cif")) - @test_throws ErrorException Crystal("SBMOF-1_overlap_diff_atom.cif", remove_duplicates=true) # duplicate but diff atom, so not repairable -end -end diff --git a/test/distance.jl b/test/distance.jl deleted file mode 100644 index 0b0051c92..000000000 --- a/test/distance.jl +++ /dev/null @@ -1,113 +0,0 @@ -module Distance_Test - -using PorousMaterials -using Test -using LinearAlgebra - -@testset "Distance Tests" begin - dxf = [-0.8, -0.4, 0.7] - nearest_image!(dxf) - @test isapprox(dxf, [0.2, -0.4, -0.3]) - - dxf = [-0.3, -0.1, -0.9] - nearest_image!(dxf) - @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] - ) - box = unit_cube() - @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])) - 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])) - - # distance (Cartesian) - 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, 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])) - 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])) - - # distance (Tests from Avogadro measurement tool) - crystal = Crystal("distance_test.cif") - @test isapprox(distance(crystal.atoms, crystal.box, 1, 3, false), 12.334, atol=0.001) # C- Ca - @test isapprox(distance(crystal.atoms, crystal.box, 1, 2, false), 14.355, atol=0.001) # C-S - @test isapprox(distance(crystal.atoms, crystal.box, 1, 2, true), 8.841, atol=0.001) # C-S - @test isapprox(distance(crystal.atoms, crystal.box, 2, 3, false), 6.292, atol=0.001) # S-Ca - - # 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) - @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) - @test o_flag - @test o_ids == [(1, 4)] - o_flag, o_ids = overlap(f, box, false, tol=0.03) - @test !o_flag - @test o_ids == [] - - # test distance function (via Avogadro) - crystal = Crystal("simple_test.cif") - @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) - @test isapprox(distance(crystal.atoms, crystal.box, 1, 5, false), 17.279, atol=0.001) - @test isapprox(distance(crystal.atoms, crystal.box, 1, 5, true), 1.531, atol=0.001) - - ### - # remove duplicates - ### - atoms = Crystal("SBMOF-1.cif").atoms - box = Crystal("SBMOF-1.cif").box - # no duplicates in SBMOF-1 - atoms_dm = remove_duplicates(atoms, box, true) - @test isapprox(atoms, atoms_dm) - # if atoms overlap but not same species, not duplicate - atoms.coords.xf[:, 5] = atoms.coords.xf[:, end] # atoms 5 and end now overlap, but 5 is C and end is Ca - atoms_dm = remove_duplicates(atoms, box, true) - @test isapprox(atoms, atoms_dm) - # finally introduce overlap between atom 5, 6, 7, all C - atoms.coords.xf[:, 6] = atoms.coords.xf[:, 5] # atoms 5 and 7 now overlap - atoms.coords.xf[:, 7] = atoms.coords.xf[:, 5] # atoms 5 and 7 now overlap - atoms_dm = remove_duplicates(atoms, box, true) - @test atoms_dm.n == atoms.n - 2 - @test isapprox(atoms[1:5], atoms_dm[1:5]) - @test isapprox(atoms[8:end], atoms_dm[6:end]) - - ### - # pairwise distances - ### - coords = Crystal("distance_tester.cif").atoms.coords - box = Crystal("distance_tester.cif").box - pd = PorousMaterials.pairwise_distances(coords, box, true) - @test isapprox(pd[2, 2], 0.0) # zero on diag - @test isapprox(pd[1, 2], pd[2, 1]) # symmetry - @test isapprox(pd[1, 2], 2.835, atol=0.01) # calculated in avogadro - @test isapprox(pd[3, 2], 6.031, atol=0.01) # calculated in avogadro - pd = PorousMaterials.pairwise_distances(coords, box, false) - @test isapprox(pd[3, 2], 17.31, atol=0.01) # calculated in avogadro -end -end diff --git a/test/gcmc_long.jl b/test/gcmc_long.jl index ab8ea9df4..d155b6210 100644 --- a/test/gcmc_long.jl +++ b/test/gcmc_long.jl @@ -8,10 +8,10 @@ using Printf using Statistics tests_to_run = Dict( - "ideal_gas" => false, + "ideal_gas" => true, "Xe in SBMOF-1" => true, - "CO2 in ZIF-71" => false, - "density grid" => false + "CO2 in ZIF-71" => false, # will not work in current version + "density grid" => true ) @testset "GCMC (long) tests" begin @@ -35,7 +35,6 @@ tests_to_run = Dict( for i = 1:length(fugacity) results, molecules = μVT_sim(empty_space, ideal_gas, temperature, fugacity[i], forcefield, n_burn_cycles=100000, n_sample_cycles=100000, verbose=false) - # n_burn_cycles=100000, n_sample_cycles=100000) @printf("fugacity %f, n_ig = %f, n_sim = %f\n", fugacity[i], n_ig[i], results["⟨N⟩ (molecules/unit cell)"]) @test isapprox(results["⟨N⟩ (molecules/unit cell)"], n_ig[i], rtol=0.05) end diff --git a/test/matter.jl b/test/matter.jl deleted file mode 100644 index 267e82c6f..000000000 --- a/test/matter.jl +++ /dev/null @@ -1,167 +0,0 @@ -module Matter_Test - -using PorousMaterials -using Test - -@testset "Matter Tests" begin - 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] - ) - @test length(c1) == 2 - c2 = Cart([7.0 10.0; - 8.0 11.0; - 9.0 12.0] - ) - - s1 = [:a, :b] - s2 = [:c, :d] - - ## Atoms - a1_f = Atoms(s1, f1) - a2_f = Atoms(s2, f2) - - 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)) - - a3_c = a1_c + a2_c - @test a3_c.n == a1_c.n + a2_c.n - @test a3_c.species[1:2] == s1 - @test a3_c.species[3:4] == s2 - @test a3_c.coords.x[:, 1:2] == c1.x - @test a3_c.coords.x[:, 3:4] == c2.x - - a3_f = a1_f + a2_f - @test a3_f.n == a1_f.n + a2_f.n - @test a3_f.species[1:2] == s1 - @test a3_f.species[3:4] == s2 - @test a3_f.coords.xf[:, 1:2] == f1.xf - @test a3_f.coords.xf[:, 3:4] == f2.xf - - ## Charges - q_vals_1 = [1.0, 4.0] - q_vals_2 = [0.4, 0.3] - - q1_c = Charges(q_vals_1, c1) - q2_c = Charges(q_vals_2, c2) - @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)) - - q3_c = q1_c + q2_c - @test q3_c.n == q1_c.n + q2_c.n - @test q3_c.q[1:2] == q_vals_1 - @test q3_c.q[3:4] == q_vals_2 - @test q3_c.coords.x[:, 1:2] == c1.x - @test q3_c.coords.x[:, 3:4] == c2.x - - q3_f = q1_f + q2_f - @test q3_f.n == q1_f.n + q2_f.n - @test q3_f.q[1:2] == q_vals_1 - @test q3_f.q[3:4] == q_vals_2 - @test q3_f.coords.xf[:, 1:2] == f1.xf - @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] - ) - wrap!(f) - @test isapprox(f, fw) - - q = Charges([-0.1, 0.2, -0.1], Frac(zeros(3, 3))) - @test neutral(q) - @test isapprox(net_charge(q), 0.0) - q.q[2] = 0.0 - @test ! neutral(q) - @test isapprox(net_charge(q), -0.2) - q = Charges{Frac}(0) - @test net_charge(q) == 0.0 - q = Charges{Cart}(0) - @test net_charge(q) == 0.0 - - # safe constructors - atoms = Atoms{Frac}(10) - @test all(isnan.(atoms.coords.xf)) - @test all(atoms.species .== :_) - atoms = Atoms{Cart}(10) - @test all(isnan.(atoms.coords.x)) - @test all(atoms.species .== :_) - - charges = Charges{Frac}(10) - @test all(isnan.(charges.coords.xf)) - @test all(isnan.(charges.q)) - charges = Charges{Cart}(10) - @test all(isnan.(charges.coords.x)) - @test all(isnan.(charges.q)) - - # getindex (array indexing!) - atoms = Crystal("SBMOF-1.cif").atoms - @test isapprox(atoms[3:4].coords, atoms.coords[3:4]) - @test isapprox(atoms[3:4].coords.xf, atoms.coords.xf[:, 3:4]) - @test atoms[3:4].species == atoms.species[3:4] - @test length(atoms[6:10].species) == 5 == atoms[6:10].n - - charges = Crystal("ATIBOU01_clean.cif").charges - @test isapprox(charges[3:4].coords, charges.coords[3:4]) - @test isapprox(charges[3:4].coords.xf, charges.coords.xf[:, 3:4]) - @test isapprox(charges[3:4].q, charges.q[3:4]) - @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] - ) - 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] - ) - 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] - ) - ) - -end -end diff --git a/test/misc.jl b/test/misc.jl index 150efe0f4..350ebf351 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -8,9 +8,8 @@ using Optim using Random @testset "Misc Tests" begin - am = read_atomic_masses() - @test isapprox(am[:H], 1.00794, atol=0.001) - @test isapprox(am[:Co], 58.9332, atol=0.001) + @test isapprox(rc[:atomic_masses][:H], 1.00794, atol=0.001) + @test isapprox(rc[:atomic_masses][:Co], 58.9332, atol=0.001) test_xyz_filename = "atoms_test" c = Cart([1.0 4.0; @@ -24,6 +23,6 @@ using Random @test isapprox(atoms, atoms_read) rm(test_xyz_filename * ".xyz") - @test read_cpk_colors()[:Li] == (204,128,255) + @test rc[:cpk_colors][:Li] == (204,128,255) end end diff --git a/test/molecule.jl b/test/molecule.jl index 57e37d084..fad9919da 100644 --- a/test/molecule.jl +++ b/test/molecule.jl @@ -49,7 +49,6 @@ pairwise_distances(m::Molecule{Frac}, box::Box) = [pairwise_distances(m.atoms.co molecule = Molecule("CO2") @test needs_rotations(molecule) @test has_charges(molecule) - atomic_masses = read_atomic_masses() @test molecule.species == :CO2 @test molecule.atoms.n == 3 @test molecule.atoms.species[1] == :C_CO2 diff --git a/test/paths.jl b/test/paths.jl deleted file mode 100644 index b03632ced..000000000 --- a/test/paths.jl +++ /dev/null @@ -1,23 +0,0 @@ -module Path_Test - -using PorousMaterials -using Test - -@testset "Path Tests" begin - # recommended way to change path to files - @eval PorousMaterials PATH_TO_CRYSTALS = joinpath(pwd(), "other_data", "other_crystals") - @eval PorousMaterials PATH_TO_MOLECULES = joinpath(pwd(), "other_data", "other_molecules") - @eval PorousMaterials PATH_TO_FORCEFIELDS = joinpath(pwd(), "other_data", "other_forcefields") - @eval PorousMaterials PATH_TO_GRIDS = joinpath(pwd(), "my_grids") - Crystal("other_SBMOF-1.cif") - Molecule("other_Xe") - LJForceField("other_Dreiding") - if ! isdir("my_grids") - mkdir("my_grids") - end - cp("data/grids/test_grid.cube", "my_grids/other_test_grid.cube", force=true) - read_cube("other_test_grid.cube") - - @test true # if made it this far :) -end -end diff --git a/test/runtests.jl b/test/runtests.jl index 6e2b49aa9..bc4b01a26 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,30 +1,25 @@ -# Details from http://www.stochasticlifestyle.com/finalizing-julia-package-documentation-testing-coverage-publishing/ +testfiles = [ + "molecule.jl", + "gcmc_quick.jl", + "grid.jl", + "misc.jl", + "forcefield.jl", + "vdw_energetics.jl", + "energy_utilities.jl", + "electrostatics.jl", + "mc_helpers.jl", + "eos.jl" + ] -function runtest(testfile::String) - @info "Testing $(testfile)" - try - include(testfile) - catch exception - @error "Exception in $(testfile)" exception - end -end +using Test, FIGlet + +font_num = 57 +FIGlet.render("Porous", FIGlet.availablefonts()[font_num]) +FIGlet.render("Materials", FIGlet.availablefonts()[font_num]) -testfiles = ["box.jl", - "matter.jl", - "crystal.jl", - "distance.jl", - "misc.jl", - "forcefield.jl", - "molecule.jl", - "vdw_energetics.jl", - "energy_utilities.jl", - "electrostatics.jl", - "mc_helpers.jl", - "eos.jl", - "assert_p1_symmetry.jl", - "grid.jl", - "gcmc_quick.jl", - "paths.jl" # must be last so as to not interfere with reading in files - ] +for testfile ∈ testfiles + @info "Running test/$testfile" + include(testfile) +end -[runtest(testfile) for testfile in testfiles] +@info "Tests complete!"