-
Notifications
You must be signed in to change notification settings - Fork 1
/
EquationOfStateFitting.jl
120 lines (106 loc) · 4.34 KB
/
EquationOfStateFitting.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
"""
# module EquationOfStateFitting
# Examples
```jldoctest
julia>
```
"""
module EquationOfStateFitting
using LinearAlgebra: det
using Compat: isnothing
using EquationsOfState
using EquationsOfState.Collections: EquationOfState
using EquationsOfState.NonlinearFitting: lsqfit
using EquationsOfState.Find: findvolume
using Kaleido: @batchlens
using MLStyle: @match
using QuantumESPRESSOBase: to_qe
using QuantumESPRESSOBase.Inputs.PWscf: PWscfInput, autofill_cell_parameters
using QuantumESPRESSOParsers.OutputParsers.PWscf: parse_total_energy,
parse_cell_parameters,
parse_head,
isjobdone
using Setfield: set
using Unitful: AbstractQuantity, ustrip, @u_str
using UnitfulAtomic
import ..Step
using ..SelfConsistentField: write_metadata
export update_alat_press, prepare, finish
function update_alat_press(
template::PWscfInput,
eos::EquationOfState,
pressure::Union{Real,AbstractQuantity},
)
volume = findvolume(PressureForm(), eos, pressure, (eps() * u"bohr^3", eos.v0 * 1.3))
alat = cbrt(ustrip(volume) / det(template.cell_parameters.data))
lenses = @batchlens(begin
_.system.celldm ∘ _[$1] # Get the `template`'s `system.celldm[1]` value
_.cell.press # Get the `template`'s `cell.press` value
end)
# Set the `template`'s `system.celldm[1]` and `cell.press` values with `alat` and `pressure`
return set(template, lenses, (alat, ustrip(pressure)))
end # function update_alat_press
# This is a helper function and should not be exported.
function _preset(step::Step{N}, template::PWscfInput) where {N}
lenses = @batchlens(begin
_.control.calculation # Get the `template`'s `control.calculation` value
_.control.verbosity # Get the `template`'s `control.verbosity` value
_.control.tstress # Get the `template`'s `control.tstress` value
_.control.tprnfor # Get the `template`'s `control.tprnfor` value
end)
# Set the `template`'s values with...
template = set(template, lenses, (N == 1 ? "scf" : "vc-relax", "high", true, true))
return isnothing(template.cell_parameters) ? autofill_cell_parameters(template) : template
end # function _preset
function prepare(
step::Step,
inputs::AbstractVector{<:AbstractString},
template::PWscfInput,
trial_eos::EquationOfState,
pressures::AbstractVector{<:Union{Real, AbstractQuantity}},
metadatafiles::AbstractVector{<:AbstractString} = map(x -> splitext(x)[1] * ".json", inputs),
verbose::Bool = false
)
# Check parameters
@assert(
length(inputs) == length(pressures) == length(metadatafiles),
"The inputs, pressures and the metadata files must be the same size!"
)
template = _preset(step, template)
# Write input and metadata
for (input, pressure, metadatafile) in zip(inputs, pressures, metadatafiles)
# Get a new `object` from the `template`, with its `alat` and `pressure` changed
object = update_alat_press(template, trial_eos, pressure)
write(input, to_qe(object, verbose = verbose)) # Write the `object` to a Quantum ESPRESSO input file
write_metadata(metadatafile, input, template)
end
return
end # function prepare
function finish(
::Step{N},
outputs::AbstractVector{<:AbstractString},
trial_eos::EquationOfState,
) where {N}
energies, volumes = zeros(length(outputs)), zeros(length(outputs))
for (i, output) in enumerate(outputs)
open(output, "r") do io
s = read(io, String)
isjobdone(s) || @warn("Job is not finished!")
energies[i] = parse_total_energy(s)[end]
volumes[i] = @match N begin
1 => parse_head(s)["unit-cell volume"]
2 => parse_cell_parameters(s)[end] |> det
_ => error("The step $N must be `1` or `2`!")
end
end
end
return lsqfit(EnergyForm(), trial_eos, volumes, energies)
end # function finish
# This is what the web service does
# function workflow(args)
# prepare(Step(1), inputs, template, trial_eos, pressures, metadatafiles)
# trial_eos = finish(outputs, trial_eos)
# prepare(Step(2), inputs, template, trial_eos, pressures, metadatafiles)
# return finish(outputs, trial_eos)
# end # function workflow
end