Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix energy grid w/ rotation & charges #239

Closed
wants to merge 6 commits into from
Closed

fix energy grid w/ rotation & charges #239

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
00d9066
update property name (bug fix)
eahenle Jun 24, 2024
2017cb3
update docs
eahenle Jun 24, 2024
b4f55b5
testset organization
eahenle Jun 24, 2024
628a03d
energy grid w/ charges & rotations
eahenle Jun 24, 2024
1b59add
adjust tests
eahenle Jun 24, 2024
abc161b
fix issue with inconsistent stack trace in doctesting
eahenle Jun 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 4 additions & 15 deletions docs/src/crystal.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,21 +61,10 @@ For many simulations, one needs to replicate the unit cell multiple times to cre
```jldoctest crystal
super_xtal = replicate(xtal, (2,2,2)) # Replicates the original unit cell once in each dimension
# output
Name: SBMOF-1.cif
Bravais unit cell of a crystal.
Unit cell angles α = 90.000000 deg. β = 100.897000 deg. γ = 90.000000 deg.
Unit cell dimensions a = 23.238600 Å. b = 11.133400 Å, c = 45.862400 Å
Volume of unit cell: 11651.776815 ų

# atoms = 960
# charges = 960
chemical formula: Ca₃₂C₄₄₈H₂₅₆O₁₉₂S₃₂
space Group: P1
symmetry Operations:
'x, y, z'
bonding graph:
# vertices = 960
# edges = 0
Crystal(Ca₃₂C₄₄₈H₂₅₆O₁₉₂S₃₂, periodic = TTT):
bounding_box : [ 23.2386 0 0;
6.81724e-16 11.1334 0;
-8.67001 3.33915e-15 45.0354]u"Å"
```

## finding other properties
Expand Down
9 changes: 7 additions & 2 deletions docs/src/matter.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,9 +57,14 @@ length(coords) # number of particles, (5)
`Coords` are immutable:

```jldoctest matter
coords.x = rand(3, 5)
try
coords.x = rand(3, 5)
catch exception
@error "Immutable!"
end
# output
ERROR: setfield!: immutable struct of type Cart cannot be changed
┌ Error: Immutable!
└ @ Main none:4
```
but we can manipulate the values of `Array{Float64, 2}` where coordinates (through `coords.x` or `coords.xf`) are stored:

Expand Down
4 changes: 2 additions & 2 deletions src/grid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -339,12 +339,12 @@ function energy_grid(crystal::Crystal, molecule::Molecule{Cart}, ljforcefield::L
else
boltzmann_factor_sum = 0.0
for _ ∈ 1:n_rotations
random_rotation!(molecule)
random_rotation!(molecule, crystal.box)

energy = PotentialEnergy(0.0, 0.0)
energy.vdw = vdw_energy(crystal, molecule, ljforcefield)
if charged_system
energy.coulomb = total(electrostatic_potential_energy(crystal, molecule,
energy.es = total(electrostatic_potential_energy(crystal, molecule,
eparams, eikr))
end
boltzmann_factor_sum += exp(-sum(energy) / temperature)
Expand Down
317 changes: 171 additions & 146 deletions test/grid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,170 +11,195 @@ using Random
CH4 = Molecule("CH4")
UFF = LJForceField("UFF")

@testset "Grid Tests" begin
# required number of pts
box = Box(1.0, 10.0, 5.0, π/2, π/2, π/2)
@test required_n_pts(box, 0.1) == (11, 101, 51)

# test read and write
grid = Grid(Box(0.7, 0.8, 0.9, 1.5, 1.6, 1.7), (3, 3, 3), rand(Float64, (3, 3, 3)),
:kJ_mol, [1., 2., 3.])
write_cube(grid, "test_grid.cube")
grid2 = read_cube("test_grid.cube")
@test isapprox(grid, grid2, atol=1e-5) # atol b/c loose precision when reading/writing to file

# nearest neighbor ID checker
n_pts = (10, 20, 30)
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.001]) == [1, 1, 1]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.999, 0.999, 0.999]) == [10, 20, 30]
idx = [0, 21, 31]
PorousMaterials._apply_pbc_to_index!(idx, n_pts)
@test idx == [10, 1, 1]
n_pts = (3, 3, 3) # so grid is [0, 0.5, 1.0]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.001]) == [1, 1, 1]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.999, 0.999, 0.999]) == [3, 3, 3]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.24]) == [1, 1, 1]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.26]) == [1, 1, 2]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.74, 0.001, 0.26]) == [2, 1, 2]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.76, 0.001, 0.26]) == [3, 1, 2]
@test PorousMaterials._arg_nearest_neighbor((4, 6, 11), [0.32, 0.42, 0.61]) == [2, 3, 7]

# accessibility grids
for zeolite in ["LTA", "SOD"]
crystal = Crystal(zeolite * ".cif")
write_xyz(crystal)
molecule = deepcopy(CH4)
forcefield = deepcopy(UFF)
grid = energy_grid(crystal, molecule, forcefield, resolution=5.)
@testset "Grid Tests" verbose=true begin
@testset "required_n_pts" verbose=true begin
box = Box(1.0, 10.0, 5.0, π/2, π/2, π/2)
@test required_n_pts(box, 0.1) == (11, 101, 51)
end

# endpoints included, ensure periodic since endpoints of grid pts included
# first cut out huge values. 1e46 == 1.00001e46
grid.data[grid.data .> 1000.0] .= 0.0
@test isapprox(grid.data[1, :, :], grid.data[end, :, :], atol=1e-7)
@test isapprox(grid.data[:, 1, :], grid.data[:, end, :], atol=1e-7)
@test isapprox(grid.data[:, :, 1], grid.data[:, :, end], atol=1e-7)
@testset "read/write" verbose=true begin
grid = Grid(Box(0.7, 0.8, 0.9, 1.5, 1.6, 1.7), (3, 3, 3), rand(Float64, (3, 3, 3)),
:kJ_mol, [1., 2., 3.])
write_cube(grid, "test_grid.cube")
grid2 = read_cube("test_grid.cube")
@test isapprox(grid, grid2, atol=1e-5) # atol b/c loose precision when reading/writing to file
end

accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=true)
@test nb_segments_blocked > 0
if zeolite == "LTA"
@test nb_segments_blocked == 8
end
@test porosity[:b4_blocking] > porosity[:after_blocking]
@testset "nearest neighbor ID checker" verbose=true begin
n_pts = (10, 20, 30)
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.001]) == [1, 1, 1]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.999, 0.999, 0.999]) == [10, 20, 30]
idx = [0, 21, 31]
PorousMaterials._apply_pbc_to_index!(idx, n_pts)
@test idx == [10, 1, 1]
n_pts = (3, 3, 3) # so grid is [0, 0.5, 1.0]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.001]) == [1, 1, 1]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.999, 0.999, 0.999]) == [3, 3, 3]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.24]) == [1, 1, 1]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.001, 0.001, 0.26]) == [1, 1, 2]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.74, 0.001, 0.26]) == [2, 1, 2]
@test PorousMaterials._arg_nearest_neighbor(n_pts, [0.76, 0.001, 0.26]) == [3, 1, 2]
@test PorousMaterials._arg_nearest_neighbor((4, 6, 11), [0.32, 0.42, 0.61]) == [2, 3, 7]
end

@test isapprox(crystal.box, accessibility_grid.box)
@testset "accessibility grids" verbose=true begin
for zeolite in ["LTA", "SOD"]
crystal = Crystal(zeolite * ".cif")
write_xyz(crystal)
molecule = deepcopy(CH4)
forcefield = deepcopy(UFF)
grid = energy_grid(crystal, molecule, forcefield, resolution=5.)

# endpoints included, ensure periodic since endpoints of grid pts included
# first cut out huge values. 1e46 == 1.00001e46
grid.data[grid.data .> 1000.0] .= 0.0
@test isapprox(grid.data[1, :, :], grid.data[end, :, :], atol=1e-7)
@test isapprox(grid.data[:, 1, :], grid.data[:, end, :], atol=1e-7)
@test isapprox(grid.data[:, :, 1], grid.data[:, :, end], atol=1e-7)

accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=true)
@test nb_segments_blocked > 0
if zeolite == "LTA"
@test nb_segments_blocked == 8
end
@test porosity[:b4_blocking] > porosity[:after_blocking]

if zeolite == "SOD"
@test all(.! accessibility_grid.data)
end
@test isapprox(crystal.box, accessibility_grid.box)

if zeolite == "SOD"
@test all(.! accessibility_grid.data)
end

# test accessibility by inserting random particles and writing to .xyz only if not accessible
nb_insertions = 100000
x = zeros(3, 0)
for i = 1:nb_insertions
xf = rand(3)
if accessible(accessibility_grid, xf)
x = hcat(x, crystal.box.f_to_c * xf)
@assert accessible(accessibility_grid, xf, (1, 1, 1))
# test accessibility by inserting random particles and writing to .xyz only if not accessible
nb_insertions = 100000
x = zeros(3, 0)
for i = 1:nb_insertions
xf = rand(3)
if accessible(accessibility_grid, xf)
x = hcat(x, crystal.box.f_to_c * xf)
@assert accessible(accessibility_grid, xf, (1, 1, 1))
else
@assert ! accessible(accessibility_grid, xf, (1, 1, 1))
end
end
if zeolite == "SOD"
# shldn't be any accessible insertions
@test length(x) == 0
else
@assert ! accessible(accessibility_grid, xf, (1, 1, 1))
xyzfilename = zeolite * "accessible_inertions.xyz"
write_xyz(Atoms([:CH4 for i = 1:size(x)[2]], Cart(x)), xyzfilename)
println("See ", xyzfilename)
end

# w./o blocking (nb = no blocking)
accessibility_grid_nb, nb_segments_blocked_nb, porosity_nb = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=false, block_inaccessible_pockets=false)
@test nb_segments_blocked_nb == 0
@test porosity_nb > porosity[:after_blocking]
@test isapprox(porosity_nb, porosity[:b4_blocking])
end
if zeolite == "SOD"
# shldn't be any accessible insertions
@test length(x) == 0
else
xyzfilename = zeolite * "accessible_inertions.xyz"
write_xyz(Atoms([:CH4 for i = 1:size(x)[2]], Cart(x)), xyzfilename)
println("See ", xyzfilename)
end
end

# w./o blocking (nb = no blocking)
accessibility_grid_nb, nb_segments_blocked_nb, porosity_nb = compute_accessibility_grid(crystal,
@testset "accessibility interpolator with replications" verbose=true begin
crystal = Crystal("LTA.cif")
molecule = deepcopy(CH4)
forcefield = deepcopy(UFF)
accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=true)

# replicate crystal and build accessibility grid that includes the other accessibility grid in a corner
repfactors = (2, 3, 1)
crystal = replicate(crystal, repfactors)
rep_accessibility_grid, rep_nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=false, block_inaccessible_pockets=false)
@test nb_segments_blocked_nb == 0
@test porosity_nb > porosity[:after_blocking]
@test isapprox(porosity_nb, porosity[:b4_blocking])
write_b4_after_grids=true)
@test all(accessibility_grid.data .== rep_accessibility_grid.data[1:7, 1:7, 1:7])
@test rep_nb_segments_blocked > 0
same_accessibility_repfactors = true
for i = 1:10000
xf = rand(3) # in (2, 3, 1) box
if ! (accessible(rep_accessibility_grid, xf) == accessible(accessibility_grid, xf, repfactors))
same_accessibility_repfactors = false
end
end
@test same_accessibility_repfactors
end

# test accessibility interpolator when there are replications
crystal = Crystal("LTA.cif")
molecule = deepcopy(CH4)
forcefield = deepcopy(UFF)
accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=true)

# replicate crystal and build accessibility grid that includes the other accessibility grid in a corner
repfactors = (2, 3, 1)
crystal = replicate(crystal, repfactors)
rep_accessibility_grid, rep_nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, resolution=2., energy_tol=0.0, verbose=false,
write_b4_after_grids=true)
@test all(accessibility_grid.data .== rep_accessibility_grid.data[1:7, 1:7, 1:7])
@test rep_nb_segments_blocked > 0
same_accessibility_repfactors = true
for i = 1:10000
xf = rand(3) # in (2, 3, 1) box
if ! (accessible(rep_accessibility_grid, xf) == accessible(accessibility_grid, xf, repfactors))
same_accessibility_repfactors = false
@testset "accessibility grid w./o pocket blocking" verbose=true begin
molecule = deepcopy(CH4)
forcefield = deepcopy(UFF)
for crystal in [Crystal("SBMOF-1.cif"), Crystal("CAXVII_clean.cif")]

# w./ blocking
accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, energy_tol=5.0, verbose=false,
write_b4_after_grids=true, energy_units=:kJ_mol)
@test nb_segments_blocked == 0

# w./o blocking (nb = no blocking)
accessibility_grid_nb, nb_segments_blocked_nb, porosity_nb = compute_accessibility_grid(crystal,
molecule, forcefield, energy_tol=5.0, verbose=false, energy_units=:kJ_mol,
write_b4_after_grids=true, block_inaccessible_pockets=false)

@test isapprox(porosity_nb, porosity[:b4_blocking])
@test nb_segments_blocked_nb == 0
@test isapprox(accessibility_grid, accessibility_grid_nb)
end
end
@test same_accessibility_repfactors

# SBMOF-1, CAXVILL_clean hv no pockets blocked. test accessibility grid w./o pocket blocking
molecule = deepcopy(CH4)
forcefield = deepcopy(UFF)
for crystal in [Crystal("SBMOF-1.cif"), Crystal("CAXVII_clean.cif")]
@testset "xf_to_id" verbose=true begin
n_pts = (4, 4, 4) # testing a grid with 4x4x4 voxels
@test all(xf_to_id(n_pts, [0.0001, 0.0001, 0.0001]) .== 1)
@test all(xf_to_id(n_pts, [0.9999, 0.9999, 0.9999]) .== n_pts[end])
@test all(xf_to_id(n_pts, [0.0001, 0.2400, 0.2600]) .== [1, 1, 2])
@test all(xf_to_id(n_pts, [-0.0001, 0.2400, 1.01]) .== [4, 1, 1]) # PBC check
end

@testset "update_density!" verbose=true begin
n_pts = (4, 4, 4) # testing a grid with 4x4x4 voxels
unit_box = unit_cube()
density_grid_c_co2 = Grid(unit_box, n_pts, zeros(n_pts...), :invers_A3, [0.0, 0.0, 0.0])
density_grid_o_co2 = Grid(unit_box, n_pts, zeros(n_pts...), :invers_A3, [0.0, 0.0, 0.0])
molecule = Molecule("CO2")
molecule = Frac(molecule, unit_box)
translate_to!(molecule, Frac([0.24, 0.26, 0.99]))
update_density!(density_grid_c_co2, [molecule], :C_CO2) # only updates for a single atom
@test isapprox(sum(density_grid_c_co2.data), 1.0)
@test isapprox(density_grid_c_co2.data[1, 2, 4], 1.0)

update_density!(density_grid_o_co2, [molecule], :O_CO2) # only updates for a single atom
@test isapprox(sum(density_grid_o_co2.data), 2.0)

translate_to!(molecule, Frac([-0.26, 0.26, 1.1])) # test PBCs
update_density!(density_grid_c_co2, [molecule], :C_CO2) # only updates for a single atom
@test isapprox(sum(density_grid_c_co2.data), 2.0)
@test isapprox(density_grid_c_co2.data[3, 2, 1], 1.0)
end

# w./ blocking
accessibility_grid, nb_segments_blocked, porosity = compute_accessibility_grid(crystal,
molecule, forcefield, energy_tol=5.0, verbose=false,
write_b4_after_grids=true, energy_units=:kJ_mol)
@test nb_segments_blocked == 0

# w./o blocking (nb = no blocking)
accessibility_grid_nb, nb_segments_blocked_nb, porosity_nb = compute_accessibility_grid(crystal,
molecule, forcefield, energy_tol=5.0, verbose=false, energy_units=:kJ_mol,
write_b4_after_grids=true, block_inaccessible_pockets=false)

@test isapprox(porosity_nb, porosity[:b4_blocking])
@test nb_segments_blocked_nb == 0
@test isapprox(accessibility_grid, accessibility_grid_nb)
@testset "xf to id" verbose=true begin
@test isapprox(id_to_xf((1, 1, 1), (10, 12, 14)), zeros(3))
@test isapprox(id_to_xf((10, 12, 14), (10, 12, 14)), ones(3))
@test isapprox(id_to_xf((2, 2, 2), (3, 3, 3)), [0.5, 0.5, 0.5])
end

# test xf_to_id
n_pts = (4, 4, 4) # testing a grid with 4x4x4 voxels
@test all(xf_to_id(n_pts, [0.0001, 0.0001, 0.0001]) .== 1)
@test all(xf_to_id(n_pts, [0.9999, 0.9999, 0.9999]) .== n_pts[end])
@test all(xf_to_id(n_pts, [0.0001, 0.2400, 0.2600]) .== [1, 1, 2])
@test all(xf_to_id(n_pts, [-0.0001, 0.2400, 1.01]) .== [4, 1, 1]) # PBC check

# test update_density!
unit_box = unit_cube()
density_grid_c_co2 = Grid(unit_box, n_pts, zeros(n_pts...), :invers_A3, [0.0, 0.0, 0.0])
density_grid_o_co2 = Grid(unit_box, n_pts, zeros(n_pts...), :invers_A3, [0.0, 0.0, 0.0])
molecule = Molecule("CO2")
molecule = Frac(molecule, unit_box)
translate_to!(molecule, Frac([0.24, 0.26, 0.99]))
update_density!(density_grid_c_co2, [molecule], :C_CO2) # only updates for a single atom
@test isapprox(sum(density_grid_c_co2.data), 1.0)
@test isapprox(density_grid_c_co2.data[1, 2, 4], 1.0)

update_density!(density_grid_o_co2, [molecule], :O_CO2) # only updates for a single atom
@test isapprox(sum(density_grid_o_co2.data), 2.0)

translate_to!(molecule, Frac([-0.26, 0.26, 1.1])) # test PBCs
update_density!(density_grid_c_co2, [molecule], :C_CO2) # only updates for a single atom
@test isapprox(sum(density_grid_c_co2.data), 2.0)
@test isapprox(density_grid_c_co2.data[3, 2, 1], 1.0)

# xf to id
@test isapprox(id_to_xf((1, 1, 1), (10, 12, 14)), zeros(3))
@test isapprox(id_to_xf((10, 12, 14), (10, 12, 14)), ones(3))
@test isapprox(id_to_xf((2, 2, 2), (3, 3, 3)), [0.5, 0.5, 0.5])
@testset "charges and rotations" verbose=true begin
xtal = Crystal("zif71_bogus_charges.cif")
strip_numbers_from_atom_labels!(xtal)
molecule = Molecule("CO2")
ljforcefield = deepcopy(UFF)
grid = energy_grid(xtal, molecule, ljforcefield, resolution=3., units=:kJ_mol, temperature=300.0, n_rotations=2, verbose=false)
@test isapprox(grid.box.a, 28.5539)
@test isapprox(grid.box.α, 1.5707963267948966)
@test grid.n_pts == (11, 11, 11)
@test all(grid.n_pts .== size(grid.data))
@test isa(grid.data, Array{Float64, 3})
@test !any(isnan, grid.data)
@test grid.units == :kJ_mol
end

end
end
Loading