-
Notifications
You must be signed in to change notification settings - Fork 22
/
atomsbase.jl
164 lines (146 loc) · 6.4 KB
/
atomsbase.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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
import AtomsBase
using AtomsBase: AbstractSystem, FlexibleSystem
using Unitful
using UnitfulAtomic
import PeriodicTable
"""
Construct an AtomsBase `FlexibleSystem` from a Chemfiles Frame.
"""
function AtomsBase.FlexibleSystem(frame::Chemfiles.Frame)
atoms = map(1:length(frame), frame) do i, atom
pos = Chemfiles.positions(frame)[:, i]u"Å"
if Chemfiles.has_velocities(frame)
velocity = Chemfiles.velocities(frame)[:, i]u"Å/ps"
else
velocity = zeros(3)u"Å/ps"
end
# Collect atomic properties
atprops = Dict(
:atomic_mass => Chemfiles.mass(atom)u"u",
:atomic_symbol => Symbol(Chemfiles.name(atom)),
:atomic_number => Chemfiles.atomic_number(atom),
:charge => Chemfiles.charge(atom)u"e_au",
:covalent_radius => Chemfiles.covalent_radius(atom)u"Å",
:vdw_radius => Chemfiles.vdw_radius(atom)*u"Å",
)
for prop in Chemfiles.list_properties(atom)
symbol = Symbol(prop)
if !hasfield(AtomsBase.Atom, symbol) && !(symbol in keys(atprops))
atprops[symbol] = Chemfiles.property(atom, prop)
end
end
AtomsBase.Atom(Chemfiles.atomic_number(atom), pos, velocity; atprops...)
end
# Collect system properties
sysprops = Dict{Symbol,Any}()
for prop in Chemfiles.list_properties(frame)
if hasfield(FlexibleSystem, Symbol(prop))
continue
elseif prop == "charge"
value = Chemfiles.property(frame, "charge")
value isa AbstractString && (value = parse(Float64, value)) # Work around a bug
sysprops[:charge] = Float64(value)u"e_au"
elseif prop == "multiplicity"
value = Chemfiles.property(frame, "multiplicity")
value isa AbstractString && (value = parse(Float64, value)) # Work around a bug
sysprops[:multiplicity] = Int(value)
else
sysprops[Symbol(prop)] = Chemfiles.property(frame, prop)
end
end
# Construct flexible system
cell_shape = Chemfiles.shape(Chemfiles.UnitCell(frame))
if cell_shape == Chemfiles.Infinite
AtomsBase.isolated_system(atoms; sysprops...)
else
@assert cell_shape in (Chemfiles.Triclinic, Chemfiles.Orthorhombic)
box = collect(eachrow(Chemfiles.matrix(Chemfiles.UnitCell(frame))))u"Å"
AtomsBase.periodic_system(atoms, box; sysprops...)
end
end
"""
Construct an AtomsBase `AbstractSystem` from a Chemfiles Frame.
"""
AtomsBase.AbstractSystem(frame::Chemfiles.Frame) = FlexibleSystem(frame)
"""
Convert a Chemfiles Frame to an AtomsBase `AbstractSystem`
"""
Base.convert(::Type{AbstractSystem}, frame::Chemfiles.Frame) = AbstractSystem(frame)
"""
Convert a Chemfiles Frame to an AtomsBase `FlexibleSystem`
"""
Base.convert(::Type{FlexibleSystem}, frame::Chemfiles.Frame) = FlexibleSystem(frame)
function Base.convert(::Type{Frame}, system::AbstractSystem{D}) where {D}
D != 3 && @warn "1D and 2D systems not yet fully supported."
frame = Chemfiles.Frame()
# Cell and boundary conditions
if AtomsBase.bounding_box(system) == AtomsBase.infinite_box(D) # System is infinite
cell = Chemfiles.UnitCell(zeros(3, 3))
Chemfiles.set_shape!(cell, Chemfiles.Infinite)
Chemfiles.set_cell!(frame, cell)
else
if any(!isequal(AtomsBase.Periodic()), AtomsBase.boundary_conditions(system))
@warn("Ignoring specified boundary conditions: Chemfiles only supports " *
"infinite or all-periodic boundary conditions.")
end
box = zeros(3, 3)
for i = 1:D
box[i, 1:D] = ustrip.(u"Å", AtomsBase.bounding_box(system)[i])
end
cell = Chemfiles.UnitCell(box)
Chemfiles.set_cell!(frame, cell)
end
if any(atom -> !ismissing(AtomsBase.velocity(atom)), system)
Chemfiles.add_velocities!(frame)
end
for atom in system
# We are using the atomic_number here, since in AtomsBase the atomic_symbol
# can be more elaborate (e.g. D instead of H or "¹⁸O" instead of just "O").
# In Chemfiles this is solved using the "name" of an atom ... to which we
# map the AtomsBase.atomic_symbol.
cf_atom = Chemfiles.Atom(PeriodicTable.elements[atomic_number(atom)].symbol)
Chemfiles.set_name!(cf_atom, string(AtomsBase.atomic_symbol(atom)))
Chemfiles.set_mass!(cf_atom, ustrip(u"u", AtomsBase.atomic_mass(atom)))
@assert Chemfiles.atomic_number(cf_atom) == AtomsBase.atomic_number(atom)
for (k, v) in pairs(atom)
if k in (:atomic_symbol, :atomic_number, :atomic_mass, :velocity, :position)
continue # Dealt with separately
elseif k == :charge
Chemfiles.set_charge!(cf_atom, ustrip(u"e_au", v))
elseif k == :vdw_radius
if v != Chemfiles.vdw_radius(cf_atom)u"Å"
@warn "Atom vdw_radius in Chemfiles cannot be mutated"
end
elseif k == :covalent_radius
if v != Chemfiles.covalent_radius(cf_atom)u"Å"
@warn "Atom covalent_radius in Chemfiles cannot be mutated"
end
elseif v isa Union{Bool, Float64, String, Vector{Float64}}
Chemfiles.set_property!(cf_atom, string(k), v)
else
@warn "Ignoring unsupported property type $(typeof(v)) in Chemfiles for key $k"
end
end
pos = convert(Vector{Float64}, ustrip.(u"Å", position(atom)))
if ismissing(AtomsBase.velocity(atom))
Chemfiles.add_atom!(frame, cf_atom, pos)
else
vel = convert(Vector{Float64}, ustrip.(u"Å/ps", AtomsBase.velocity(atom)))
Chemfiles.add_atom!(frame, cf_atom, pos, vel)
end
end
for (k, v) in pairs(system)
if k in (:bounding_box, :boundary_conditions)
continue # Already dealt with
elseif k in (:charge, )
Chemfiles.set_property!(frame, string(k), Float64(ustrip(u"e_au", v)))
elseif k in (:multiplicity, )
Chemfiles.set_property!(frame, string(k), Float64(v))
elseif v isa Union{Bool, Float64, String, Vector{Float64}}
Chemfiles.set_property!(frame, string(k), v)
else
@warn "Ignoring unsupported property type $(typeof(v)) in Chemfiles for key $k"
end
end
frame
end