Add Hessian calculation support

Project.toml

@@ -1,7 +1,7 @@
 name = "AtomsCalculators"
 uuid = "a3e0e189-c65a-42c1-833c-339540406eb1"
 authors = ["JuliaMolSim contributors"]
-version = "0.1.1"
+version = "0.1.2-DEV"
 
 [deps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"

docs/src/interface-definition.md

@@ -1,26 +1,28 @@
 # Interface Definition
 
-There are two alternative ways to call the interface: using functions `potential_energy`, `forces` and [virial](https://en.wikipedia.org/wiki/Virial_stress), or using `calculate`
-function together with `Energy`, `Forces` and `Virial`.
+There are two alternative ways to call the interface: using functions `potential_energy`, `forces`, [virial](https://en.wikipedia.org/wiki/Virial_stress) and `hessian`, or using `calculate`
+function together with `Energy`, `Forces`, `Virial` and `Hessian`.
 
 Individual calls are implemented by dispatching `AtomsCalculators` functions
 
 - `AtomsCalculators.potential_energy` for potential energy calculation
 - `AtomsCalculators.forces` for allocating force calculation and/or...
 - `AtomsCalculators.forces!` for non-allocating force calculation
 - `AtomsCalculators.virial` for [virial](https://en.wikipedia.org/wiki/Virial_stress) calculation
+- `AtomsCalculators.hessian` for Hessian
 
 The `calculate` interface is implemented by dispatching to
 
 - `AtomsCalculators.calculate` using `AtomsCalculators.Energy()` as the first argument for energy calculation
 - `AtomsCalculators.calculate` using `AtomsCalculators.Forces()` as the first argument for forces calculaton
 - `AtomsCalculators.calculate` using `AtomsCalculators.Virial()` as the first argument for virial calculation
+- `AtomsCalculators.calculate` using `AtomsCalculators.Hessian()` as the first argument for Hessian calculation
 
 You do not need to implement all of these by yourself. There is macro that will help implement the other calls. 
 
 Each of the individual calls have two common inputs: `AtomsBase.AbstractSystem` compatible structure and a `calculator` that incudes details of the calculation method. Calculate interface has additionally the type of calculation as the first input. You can tune calculation by passing keyword arguments, which can be ignored, but they need to be present in the function definition.
 
-`potential_energy`, `forces`, `forces!` and `virial`:
+`potential_energy`, `forces`, `forces!`, `virial` and `hessian`:
 
 - First input is `AtomsBase.AbstractSystem` compatible structure
 - Second input is `calculator` structure
@@ -29,7 +31,7 @@ Each of the individual calls have two common inputs: `AtomsBase.AbstractSystem`
 
 `calculate`:
 
-- First input is either `Energy()`, `Forces()` or `Virial()`
+- First input is either `Energy()`, `Forces()`, `Virial()` or `Hessian()`
 - Second is `AtomsBase.AbstractSystem` compatible structure
 - Third is `calculator` structure
 - Method has to accept keyword arguments (they can be ignored)
@@ -41,6 +43,7 @@ Outputs for the functions need to have following properties
 - Energy is a subtype of `Number` that has a unit with dimensions of energy (mass * length^2 / time^2)
 - Force output is a subtype of `AbstractVector` with element type also a subtype of AbstractVector (length 3 in 3D) and unit with dimensions of force (mass * length / time^2). With additional property that it can be reinterpret as a matrix
 - Virial is a square matrix (3x3 in 3D) that has units of force times length or energy
+- Hessian is a square matrix, with a side length of dimension times number of atoms and unit of energy divided by square of distance
 - Calculate methods return a [NamedTuple](https://docs.julialang.org/en/v1/base/base/#Core.NamedTuple) that uses keys `:energy`, `:forces` and `:virial` to identify the results, which have the types defined above
 
 
@@ -158,6 +161,40 @@ AtomsCalculators.@generate_interface function AtomsCalculators.forces!(f::Abstra
 end
 ```
 
+### Implementing Hessian
+
+Hessian can be implemented in similar ways. Using either direct call
+
+```julia
+AtomsCalculators.@generate_interface function AtomsCalculators.hessian(system, calculator::MyType; kwargs...)
+    # we can ignore kwargs... or use them to tune the calculation
+    # or give extra information like pairlist
+
+    # add your own definition
+    l = length(system) * length(position(system, 1))
+    return zeros(l,l) * (u"eV" / u"Å"^2)
+end
+```
+
+or calculator call
+
+```julia
+AtomsCalculators.@generate_interface function AtomsCalculators.calculate(
+    ::AtomsCalculators.Hessian, 
+    system, 
+    calculator::MyTypeC; 
+    kwargs...
+)
+    # we can ignore kwargs... or use them to tune the calculation
+    # or give extra information like pairlist
+
+    # add your own definition
+    l = length(system) * length(position(system, 1))
+    h = zeros(l,l) * (u"eV" / u"Å"^2)
+    return ( hessian = h, )
+end
+```
+
 ## Other Automatically Generated Calls
 
 Many methods have optimized calls when energy and forces (and/or virial) are calculated together. To allow access to these calls there are also calls
@@ -188,11 +225,12 @@ end
 
 Defining this does not overload the corresponding non-allocating call - you need to do that separately.
 
-Output for the combination methods is defined to have keys `:energy`, `:forces` and `:virial`. You can access them with
+Output for the combination methods is defined to have keys `:energy`, `:forces`, `:virial` and `:hessian`. You can access them with
 
 - `output[:energy]` for energy
 - `output[:forces]` for forces
 - `output[:virial]` for viral
+- `output[:hessian]` for hessian
 
 The type of the output can be [NamedTuple](https://docs.julialang.org/en/v1/base/base/#Core.NamedTuple), as was in the example above, or any structure that has the keys implemented and also supports splatting. The reason for this is that this allows everyone to implement the performance functions without being restricted to certain output type, and to allow using `haskey` to check the output.
 
@@ -203,6 +241,7 @@ We have implemented function calls to help you testing the API. There is one cal
 - `test_potential_energy` to test potential_energy call
 - `test_forces` to test both allocating and non-allocating force calls
 - `test_virial` to test virial call
+- `test_hessian` to test hessian call
 
 To get these access to these functions you need to call
 

src/interface.jl

@@ -2,6 +2,7 @@
 struct Energy end
 struct Forces end
 struct Virial end
+struct Hessian end
 
 function potential_energy end 
 
@@ -13,6 +14,8 @@ function virial end
 
 function calculate end
 
+function hessian end
+
 promote_force_type(::Any, ::Any) = SVector(1., 1., 1.) * u"eV/Å" |> typeof
 
 zero_forces(system, calc) = zeros( promote_force_type(system, calc), length(system) )

src/submodules/AtomsCalculatorsTesting.jl

@@ -7,6 +7,7 @@ using Test
 export test_potential_energy
 export test_forces
 export test_virial
+export test_hessian
 
 
 """
@@ -119,4 +120,23 @@ end
 
 
 
+function test_hessian(sys, calculator; kwargs...)
+    @testset "Test hessian for $(typeof(calculator))" begin
+        l = length(sys) * length(position(sys, 1))
+        h = AtomsCalculators.hessian(sys, calculator; kwargs...)
+        @test isa(h, AbstractArray)
+        @test eltype(h) <: Number
+        @test size(h) == (l,l)
+        h_cpu_array = Array(h)
+        @test dimension(h_cpu_array[1,1]) == dimension(u"J/m^2")
+        h2 = AtomsCalculators.hessian(sys, calculator; dummy_kword235887=2, kwargs...)
+        @test all( h .≈ h2 )
+        hc = AtomsCalculators.calculate(AtomsCalculators.Hessian(), sys, calculator; kwargs...)
+        @test isa(hc, NamedTuple)
+        @test haskey(hc, :hessian)
+        @test all( h .≈ hc[:hessian] )
+    end
+end
+
+
 end

src/utils.jl

@@ -57,6 +57,8 @@ macro generate_interface(expr)
             q = generate_forces_from_calculator( calculator_type )
         elseif type[:type] == :Virial
             q = generate_virial( calculator_type )
+        elseif type[:type] == :Hessian
+            q = generate_hessian( calculator_type )
         end
     else
         if type[:type] == :potential_energy
@@ -77,6 +79,8 @@ macro generate_interface(expr)
             end
         elseif type[:type] == :virial
             q = generate_calculator_virial( calculator_type )
+        elseif type[:type] == :hessian
+            q = generate_calculator_hessian( calculator_type )
         end
     end
     ex = quote
@@ -115,14 +119,14 @@ function determine_type_calculation(expr)
         error("function definition is not correct")
     end
     type_of_calculation = expr.args[1].args[1].args[end].value
-    if type_of_calculation in [:potential_energy, :forces, :forces!, :virial]
+    if type_of_calculation in [:potential_energy, :forces, :forces!, :virial, :hessian]
         return (;
             :type => type_of_calculation,
             :calculator => false
         )
     elseif type_of_calculation == :calculate # calculator interface -> find calc type too
         # Need to have definition of AtomsCalculators.Energy() etc.
-        if expr.args[1].args[3].args[1].args[2].value in [:Energy, :Forces, :Virial]
+        if expr.args[1].args[3].args[1].args[2].value in [:Energy, :Forces, :Virial, :Hessian]
             return (;
                 :type => expr.args[1].args[3].args[1].args[2].value,
                 :calculator => true
@@ -229,4 +233,25 @@ function generate_forces_from_calculator(calc_type)
         $q1
         $q2
     end
+end
+
+
+function generate_calculator_hessian(calc_type)
+    q = quote 
+        function AtomsCalculators.calculate(::AtomsCalculators.Hessian, system, calculator::$calc_type; kwargs...)
+            h = AtomsCalculators.hessian(system, calculator; kwargs...)
+            return ( hessian = h, )
+        end
+    end
+    return q
+end
+
+function generate_hessian(calc_type)
+    q = quote
+        function AtomsCalculators.hessian(system, calculator::$calc_type; kwargs...)
+            h = AtomsCalculators.calculate(AtomsCalculators.Hessian(), system, calculator; kwargs...)
+            return h[:hessian]
+        end
+    end
+    return q
 end

test/runtests.jl

@@ -40,6 +40,16 @@ using AtomsCalculators.AtomsCalculatorsTesting
         # add your own definition
         return AtomsCalculators.zero_forces(system, calculator)
     end
+
+    AtomsCalculators.@generate_interface function AtomsCalculators.hessian(system, calculator::MyType; kwargs...)
+        # we can ignore kwargs... or use them to tune the calculation
+        # or give extra information like pairlist
+
+        # add your own definition
+        l = length(system) * length(position(system, 1))
+        return zeros(l,l) * (u"eV" / u"Å"^2)
+    end
+
 
     AtomsCalculators.@generate_interface function AtomsCalculators.forces!(f::AbstractVector, system, calculator::MyOtherType; kwargs...)
         @assert length(f) == length(system)
@@ -95,6 +105,21 @@ using AtomsCalculators.AtomsCalculatorsTesting
         return ( forces = f, )
     end
 
+    AtomsCalculators.@generate_interface function AtomsCalculators.calculate(
+        ::AtomsCalculators.Hessian, 
+        system, 
+        calculator::MyTypeC; 
+        kwargs...
+    )
+        # we can ignore kwargs... or use them to tune the calculation
+        # or give extra information like pairlist
+
+        # add your own definition
+        l = length(system) * length(position(system, 1))
+        h = zeros(l,l) * (u"eV" / u"Å"^2)
+        return ( hessian = h, )
+    end
+
     hydrogen = isolated_system([
     :H => [0, 0, 0.]u"Å",
     :H => [0, 0, 1.]u"Å"
@@ -103,11 +128,13 @@ using AtomsCalculators.AtomsCalculatorsTesting
     test_potential_energy(hydrogen, MyType())
     test_forces(hydrogen, MyType())
     test_virial(hydrogen, MyType())
+    test_hessian(hydrogen, MyType())
     test_forces(hydrogen, MyOtherType())
 
     test_potential_energy(hydrogen, MyTypeC())
     test_forces(hydrogen, MyTypeC())
     test_virial(hydrogen, MyTypeC())
+    test_hessian(hydrogen, MyTypeC())
 
     efv = AtomsCalculators.energy_forces_virial(hydrogen, MyType())
     @test haskey(efv, :energy)