-
Notifications
You must be signed in to change notification settings - Fork 2
/
guessbonds.jl
97 lines (91 loc) · 4.81 KB
/
guessbonds.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
86
87
88
89
90
91
92
93
94
95
96
97
# Bond guessing algorithm derived from from Chemfiles, itself from VMD
# Data from Blue Obelisk's data repository, corrected for H, C, O, N, S and F
# (see https://github.com/chemfiles/chemfiles/issues/301#issuecomment-574100048)
const vdwradii = Float32[1.0, 1.4, 2.2, 1.9, 1.8, 1.5, 1.4, 1.3, 1.2, 1.54, 2.4,
2.2, 2.1, 2.1, 1.95, 1.9, 1.8, 1.88, 2.8, 2.4, 2.3,
2.15, 2.05, 2.05, 2.05, 2.05, 2.0, 2.0, 2.0, 2.1, 2.1,
2.1, 2.05, 1.9, 1.9, 2.02, 2.9, 2.55, 2.4, 2.3, 2.15,
2.1, 2.05, 2.05, 2.0, 2.05, 2.1, 2.2, 2.2, 2.25, 2.2,
2.1, 2.1, 2.16, 3.0, 2.7, 2.5, 2.48, 2.47, 2.45, 2.43,
2.42, 2.4, 2.38, 2.37, 2.35, 2.33, 2.32, 2.3, 2.28,
2.27, 2.25, 2.2, 2.1, 2.05, 2.0, 2.0, 2.05, 2.1, 2.05,
2.2, 2.3, 2.3, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.4, 2.0,
2.3, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0,
2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
# Data from PeriodicTable.jl
const ismetal = Bool[0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1]
const ismetalloid = Bool[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
"""
guess_bonds(pos, types, mat, options)
Return the bonds guessed from the positions, types and cell matrix, given as a
`Vector{Vector{Tuple{Int,Float32}}}`.
The `i`-th entry of the list is a list, whose entries are of the form `(j, dist)` which
indicates that the representatives of vertices `i` and `j` distant of at most `dist` are
bonded together.
"""
function guess_bonds(pos, types, mat, options)
# Algorithm from chemfiles, itself from VMD
n = length(pos)
bonds = [Tuple{Int,Float32}[] for _ in 1:n]
options.bonding == Bonding.NoBond && return bonds
@ifwarn if options.bonding != Bonding.Guess
@warn "Guessing bonds with custom algorithm (from Chemfiles and VMD). This may take a while for big structures and may be inexact."
@info "To avoid guessing bonds, use a file format that contains the bonds."
end
@toggleassert n == length(types)
radii = Vector{Float32}(undef, n)
for (i, typ) in enumerate(types)
t = get(atomic_numbers, typ, nothing)
if t isa Int
radii[i] = vdwradii[t]*(1 + options.wider_metallic_bonds*(ismetal[t]|ismetalloid[t])*0.5)
else
radii[i] = 0.0
@ifwarn @warn lazy"Unrecognized atom type \"$t\" will be considered a dummy atom."
end
end
cutoff_sq = (3*(options.cutoff_coeff^3.1) * max(maximum(radii), 0.833))^2
cutoff2 = 13*options.cutoff_coeff/15
pd2 = PeriodicDistance2(mat)
for i in 1:n
radius_i = radii[i]
iszero(radius_i) && continue
posi = pos[i]
typi = types[i]
#ignoreifmetallic = typi === :C
#ignoreifC = ismetal[atomic_numbers[typi]]
skiphomoatomic = typi ∈ options.ignore_homoatomic_bonds ||
(options.ignore_homometallic_bonds && ismetal[atomic_numbers[typi]])
acceptonlyO = options.structure === StructureType.Zeolite && typi !== :O
acceptallbutO = options.structure === StructureType.Zeolite && typi === :O
for j in (i+1):n
typj = types[j]
skiphomoatomic && typj === typi && continue
acceptonlyO && typj !== :O && continue
acceptallbutO && typj === :O && continue
#(ignoreifC & (typj === :C)) && continue
#ignoreifmetallic && ismetal[atomic_numbers[typj]] && continue
radius_j = radii[j]
iszero(radius_j) && continue
posj = pos[j]
d2 = pd2(posi, posj)
maxdist = cutoff2*(radius_i + radius_j)
if d2 < cutoff_sq && 0.5^2 < d2 < maxdist^2
push!(bonds[i], (j, maxdist))
push!(bonds[j], (i, maxdist))
end
end
end
return bonds
end