-
Notifications
You must be signed in to change notification settings - Fork 11
/
energy_min.jl
85 lines (72 loc) · 3.71 KB
/
energy_min.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
"""
minimized_molecule, min_energy = find_energy_minimum(xtal, molecule, ljff) # molecule set at initial guess
find the minimum energy position, and associated minimum energy, of a molecule in a crystal.
n.b. if molecule has more than one atom, it will *not* minimize over the orientation (rotations).
the optimizer needs an initial estimate of the minimum energy position.
pass molecule with good initial position.
if you don't have a good initial position, use [`find_energy_minimum_gridsearch`](@ref).
# Arguments
- `xtal::Crystal`: the crystal
- `molecule::Molecule`: the molecule, whose position we seek to tune until we reach a local minimum. must start at a good initial position close to the minimum.
- `ljff::LJForceField`: the force field used to calculate crystal-molecule interaction energies
# 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(xtal::Crystal,
molecule::Molecule,
ljff::LJForceField,
)
if needs_rotations(molecule)
@warn "needs rotations. does not optimize over configurations, only over center of mass"
end
if has_charges(molecule) && has_charges(xtal)
error("cannot handle electrostatics")
end
# make sure replication factors sufficient
rep_factors = replication_factors(xtal, sqrt(ljff.r²_cutoff))
xtal = replicate(xtal, rep_factors)
# make sure molecule is in fractional coords
if isa(molecule, Molecule{Cart})
molecule = Frac(molecule, xtal.box)
else
translate_to!(molecule, Frac(molecule.com.xf ./ rep_factors))
end
# xf::Array{Float64, 1}
function energy(xf)
# make sure the coords are fractional
xf = mod.(xf, 1.0)
# move probe molecule
translate_to!(molecule, Frac(xf))
# calculate the guest-host VDW interaction and return
return vdw_energy(xtal, molecule, ljff) * 8.314 / 1000 # units: kJ/mol
end
# find minimum energy position
res = optimize(energy, molecule.com.xf)
xf_min = mod.(mod.(res.minimizer, 1.0) .* rep_factors, 1.0)
# translate molecule to min energy position to return it.
translate_to!(molecule, Frac(xf_min))
return molecule, res.minimum
end
"""
xf₀ = find_energy_minimum_gridsearch(xtal, molecule, ljff; n_pts=(50, 50, 50))
perform an [`energy_grid`](@ref) calculation and, via a grid search, find the minimum energy position of a molecule.
# Arguments
- `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
- `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;
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; resolution=resolution)
E_min, idx_min = findmin(grid.data)
xf_minE = id_to_xf(Tuple(idx_min), grid.n_pts) # fractional coords of min
molecule = Frac(molecule, xtal.box)
translate_to!(molecule, Frac(xf_minE))
return molecule, E_min
end