diff --git a/CHANGELOG.md b/CHANGELOG.md
index efc5d6f..d200572 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
 # CHANGES
 
+## v1.6.1 May 18, 2026
+- fixed bug in `PointEvaluator` initialization which broke automatic differentiation of the `evaluate!` function with respect to spatial coordinates
+
 ## v1.6.0 April 28, 2026
   - unicode plotting extension marked as deprecated (will be removed in next major version), they can be realized now with GridVisualize 1.9 and backend Plotter = UnicdePlots
   - added scalarplot, vectorplot and broken_scalarplot defined for FEVectorBlocks (moved from ExtendableFEM)
diff --git a/Project.toml b/Project.toml
index a09366c..7294e7a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ExtendableFEMBase"
 uuid = "12fb9182-3d4c-4424-8fd1-727a0899810c"
-version = "1.6.0"
+version = "1.6.1"
 authors = ["Christian Merdon <merdon@wias-berlin.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de>"]
 
 [deps]
@@ -52,10 +52,11 @@ julia = "1.9"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 Term = "22787eb5-b846-44ae-b979-8e399b8463ab"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [targets]
-test = ["Aqua", "ExampleJuggler", "ExplicitImports", "GridVisualize", "Term", "Test", "UnicodePlots"]
+test = ["Aqua", "ExampleJuggler", "ExplicitImports", "ForwardDiff", "GridVisualize", "Term", "Test", "UnicodePlots"]
diff --git a/src/point_evaluator.jl b/src/point_evaluator.jl
index 0715dda..2aa5743 100644
--- a/src/point_evaluator.jl
+++ b/src/point_evaluator.jl
@@ -146,7 +146,7 @@ function initialize!(O::PointEvaluator{T, TCoeff, UT}, sol; time = 0, kwargs...)
     O.L2G = []
     for EG in EGs
         ## FE basis evaluator for EG
-        push!(O.BE_args, [FEEvaluator(FES_args[j], O.ops_args[j], QuadratureRule{T, EG}(0); AT = AT) for j in 1:nargs])
+        push!(O.BE_args, [FEEvaluator(FES_args[j], O.ops_args[j], QuadratureRule{T, EG}(0); T, AT) for j in 1:nargs])
 
         ## L2G map for EG
         push!(O.L2G, L2GTransformer(EG, xgrid, gridAT))
diff --git a/test/runtests.jl b/test/runtests.jl
index b5b793c..76892f6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,8 @@ using ExtendableFEMBase
 using ExtendableSparse
 using ExplicitImports
 using ExampleJuggler
+using ForwardDiff
+using LinearAlgebra
 using SparseArrays
 using Aqua
 
@@ -143,6 +145,37 @@ function exact_function(::Val{3}, polyorder)
     return polynomial, exact_integral, gradient, hessian
 end
 
+function point_evaluator_closure(coeffs = [1.0, 0.0, 0.0])
+    input_types = Dict{DataType, Any}()
+    result_types = Dict{DataType, Any}()
+
+    grid_types = Dict{DataType, Any}()
+    FES_types = Dict{DataType, Any}()
+    RVec_types = Dict{DataType, Any}()
+    PE_types = Dict{DataType, Any}()
+
+    function closure(x::Vector{T}) where {T}
+        if !haskey(PE_types, T)
+            input_types[T] = zeros(T, 2)
+            result_types[T] = zeros(T, 1)
+
+            grid_types[T] = reference_domain(Triangle2D, T)
+            FES_types[T] = FESpace{H1P1{1}}(grid_types[T])
+            RVec_types[T] = FEVector(FES_types[T])
+            RVec_types[T].entries .= coeffs
+            PE_types[T] = PointEvaluator([(1, Identity)], [RVec_types[T][1]]; Tv = T, TCoeff = T)
+        end
+
+        input_types[T][1] = log(x[1])
+        input_types[T][2] = log(x[2])
+
+        ExtendableFEMBase.evaluate!(result_types[T], PE_types[T], input_types[T])
+
+        return result_types[T]
+    end
+
+    return closure
+end
 
 function run_all_tests()
     return begin
diff --git a/test/test_pointevaluator.jl b/test/test_pointevaluator.jl
index 8602ec7..7bac5c9 100644
--- a/test/test_pointevaluator.jl
+++ b/test/test_pointevaluator.jl
@@ -6,6 +6,7 @@ function run_pointevaluator_tests()
         println("======================")
         test_pointevaluation2D()
         test_pointevaluation3D()
+        test_pointevaluation_differentiability()
     end
 end
 
@@ -52,3 +53,25 @@ function test_pointevaluation3D()
     evaluate!(eval, PE, x)
     return @test abs(eval[1] - sum(x)) < 1.0e-15
 end
+
+function test_pointevaluation_differentiability()
+    ## check if point evaluation can be differentiated
+    ## wrt. coordinates
+    coefficient_test_cases = I(3)
+    x = [exp(1.0), exp(0.0)]
+
+    evals = [[0.0], [1.0], [0.0]]
+    derivs = [
+        [-1 / exp(1.0) -1.0;],
+        [1 / exp(1.0) 0.0;],
+        [0.0 1.0;],
+    ]
+
+    for (case_index, coeffs) in enumerate(eachcol(coefficient_test_cases))
+        f = point_evaluator_closure(coeffs)
+        @test norm(f(x) - evals[case_index], Inf) < tolerance
+        @test norm(ForwardDiff.jacobian(f, x) - derivs[case_index], Inf) < tolerance
+    end
+
+    return nothing
+end