Add rudimentary (untested) QE calculator

Project.toml

@@ -7,6 +7,10 @@ version = "0.1.0"
 ASEconvert = "3da9722f-58c2-4165-81be-b4d7253e8fd2"
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 DFTK = "acf6eb54-70d9-11e9-0013-234b7a5f5337"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+QuantumEspresso_jll = "74603b90-2fcf-5710-a7d7-830b31b8b33c"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 
 [compat]
 julia = "1.6"

src/ElectronicStructure.jl

@@ -3,7 +3,10 @@ module ElectronicStructure
 export AbstractCalculator, AbstractState, AbstractParameters
 export DftkCalculator, DftkParameters, DftkState, calculate, energy
 
-# Write your package code here.
-include("sandbox.jl")
+using AtomsBase
+
+include("interface.jl")
+include("dftk.jl")
+include("quantum_espresso.jl")
 
 end

src/sandbox.jl → src/dftk.jl

@@ -1,9 +1,4 @@
 import DFTK
-using AtomsBase
-
-abstract type AbstractCalculator end
-abstract type AbstractState end
-abstract type AbstractParameters end
 
 struct DftkCalculator <: AbstractCalculator end
 
@@ -33,7 +28,7 @@ function calculate(calc::DftkCalculator, params::DftkParameters)
     calculate(calc, DftkState(params))
 end
 
-function calculate(calc::DftkCalculator, state::DftkState)
+function calculate(::DftkCalculator, state::DftkState)
     scfres = state.scfres
     scfres = DFTK.self_consistent_field(state.basis;
                                         scfres.ψ, scfres.ρ, state.params.scf_kwargs...)

src/interface.jl

@@ -0,0 +1,7 @@
+# Crude definition of interface
+abstract type AbstractCalculator end
+abstract type AbstractState      end
+abstract type AbstractParameters end
+
+function calculate end
+function energy    end

src/quantum_espresso.jl

@@ -0,0 +1,93 @@
+using ASEconvert
+using PythonCall
+using QuantumEspresso_jll
+using MPI
+using Unitful
+using UnitfulAtomic
+
+
+Base.@kwdef struct QeCalculator <: AbstractCalculator
+    # TODO QE setup parameters
+end
+
+Base.@kwdef struct QeParameters <: AbstractParameters
+    # TODO Keywords currently based on ASE
+    system::AbstractSystem
+    ecutwfc     = 40
+    conv_thr    = 1e-11
+    tstress     = true
+    tprnfor     = true
+    smearing    = "gaussian"
+    mixing_mode = "plain"
+    mixing_beta = 0.7
+    mixing_ndim = 10
+    kpts        = (1, 1, 1)
+    occupations = "smearing"
+    degauss     = 0.01
+    input_dft   = "pbe"
+    electron_maxstep  = 100
+    pseudopotentials  = Dict{String,String}()
+    extra_parameter   = Dict{Symbol,Any}()
+    working_directory = mktempdir(pwd())
+    n_mpi_procs       = MPI.Comm_size(MPI.COMM_WORLD)
+    n_threads         = BLAS.get_num_threads()
+end
+
+
+
+function convert(::Type{QeParameters}, params::DftkParameters)
+    # Convert DFTK parameters to QE parameters
+    #
+    # keep in mind the unit conversion
+    error("TODO")
+end
+
+
+struct QeState <: AbstractState
+    params::QeParameters
+    ase_atoms::Py
+end
+
+function QeState(params::QeParameters)
+    ase_atoms = convert_ase(params.system)
+    ase_atoms.calc = pyimport("ase.calculators.espresso").Espresso(;
+        label="espresso",
+        params.input_dft,
+        params.pseudopotentials,
+        params.kpts,
+        params.ecutwfc,
+        params.tstress,
+        params.tprnfor,
+        params.mixing_mode,
+        params.mixing_beta,
+        params.conv_thr,
+        params.occupations,
+        params.smearing,
+        params.degauss,
+        params.electron_maxstep,
+        params.mixing_ndim,
+        params.extra_parameter...
+    )
+    QeState(params, ase_atoms)
+end
+
+function calculate(calc::QeCalculator, params::QeParameters)
+    calculate(calc, QeState(params))
+end
+
+function calculate(::QeCalculator, state::QeState)
+    n_mpi_procs = state.params.n_mpi_procs
+    MPI.mpiexec() do mpirun
+        QuantumEspresso_jll.pwscf() do pwscf
+            qe_command = "$mpirun -np $n_mpi_procs $pwscf -in PREFIX.pwi > PREFIX.pwo"
+            state.ase_atoms.calc.command = qe_command
+            withenv("OMP_NUM_THREADS" => state.params.n_threads) do
+                state.ase_atoms.calc.get_potential_energy()
+            end
+        end
+    end
+end
+
+function energy(state::QeState)
+    austrip(state.ase_calculator.get_potential_energy() * u"eV")
+end