diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index a3b1e51..6096e30 100644
--- a/src/LatticeGPU.jl
+++ b/src/LatticeGPU.jl
@@ -24,6 +24,7 @@ include("Space/Space.jl")
 using .Space
 export SpaceParm
 export up, dw, updw, global_point
+export BC_PERIODIC, BC_OPEN, BC_SF_AFWB, BC_SF_ORBI
 
 include("Fields/Fields.jl")
 using .Fields
diff --git a/src/Space/Space.jl b/src/Space/Space.jl
index 9c18850..012efd5 100644
--- a/src/Space/Space.jl
+++ b/src/Space/Space.jl
@@ -14,7 +14,12 @@ module Space
 
 import Base.show
 
-struct SpaceParm{N,M,D}
+const BC_PERIODIC = 0
+const BC_SF_ORBI  = 1
+const BC_SF_AFWB  = 2
+const BC_OPEN     = 3
+
+struct SpaceParm{N,M,B,D}
     ndim::Int64
     iL::NTuple{N,Int64}
     npls::Int64
@@ -51,21 +56,58 @@ struct SpaceParm{N,M,D}
         end
 
         D = prod(y)
-        return new{N,M,D}(N, x, M, tuple(pls...), y,
+        return new{N,M,BC_PERIODIC,D}(N, x, M, tuple(pls...), y,
+                        tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r))
+    end
+
+    function SpaceParm{N}(x, y, ibc) where {N}
+        M = convert(Int64, round(N*(N-1)/2))
+        N == length(x) || throw(ArgumentError("Lattice size incorrect length for dimension $N"))
+        N == length(y) || throw(ArgumentError("Block   size incorrect length for dimension $N"))
+
+        pls = Vector{Tuple{Int64, Int64}}()
+        for i in N:-1:1
+            for j in 1:i-1
+                push!(pls, (i,j))
+            end
+        end
+
+        r  = div.(x, y)
+        rS = ones(N)
+        yS = ones(N)
+        for i in 2:N
+            for j in 1:i-1
+                rS[i] = rS[i]*r[j]
+                yS[i] = yS[i]*y[j]
+            end
+        end
+
+        D = prod(y)
+        return new{N,M,ibc,D}(N, x, M, tuple(pls...), y,
                         tuple(yS...), tuple(r...), tuple(rS...), prod(y), prod(r))
     end
 end
 export SpaceParm
-function Base.show(io::IO, lp::SpaceParm)
-    println(io, "Lattice dimensions: ", lp.ndim)
+function Base.show(io::IO, lp::SpaceParm{N,M,B,D}) where {N,M,B,D}
+    println(io, "Lattice dimensions:       ", lp.ndim)
 
-    print(io, "Lattice size:       ")
+    print(io, "Lattice size:             ")
     for i in 1:lp.ndim-1
         print(io, lp.iL[i], " x ")
     end
     println(io,lp.iL[end])
+    print(io, "Time boundary conditions: ")
+    if B == BC_PERIODIC
+        println(io, "PERIODIC")
+    elseif B == BC_OPEN
+        println(io, "OPEN")
+    elseif B == BC_SF_AFWB
+        println(io, "SF (AFW option B)")
+    elseif B == BC_SF_ORBI
+        println(io, "SF (orbifold improvement)")
+    end
     
-    print(io, "Thread block size:  ")
+    print(io, "Thread block size:        ")
     for i in 1:lp.ndim-1
         print(io, lp.blk[i], " x ")
     end
@@ -239,5 +281,6 @@ end
 end
     
 export up, dw, updw, global_point, point_index, point_coord, point_time
+export BC_PERIODIC, BC_OPEN, BC_SF_AFWB, BC_SF_ORBI
 
 end
diff --git a/src/YM/YMact.jl b/src/YM/YMact.jl
index bab4918..2907f27 100644
--- a/src/YM/YMact.jl
+++ b/src/YM/YMact.jl
@@ -9,7 +9,7 @@
 ### created: Mon Jul 12 18:31:19 2021
 ###                               
 
-function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,D}) where {T,N,M,D}
+function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -95,7 +95,7 @@ function krnl_impr!(plx, U::AbstractArray{T}, c0, c1, lp::SpaceParm{N,M,D}) wher
     return nothing
 end
 
-function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M,D}
+function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
     
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -132,7 +132,7 @@ function krnl_plaq!(plx, U::AbstractArray{T}, lp::SpaceParm{N,M,D}) where {T,N,M
     return nothing
 end
 
-function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceParm{N,M,D}) where {T,N,M,D}
+function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -172,7 +172,7 @@ function krnl_force_wilson_pln!(frc1, frc2, U::AbstractArray{T}, ipl, lp::SpaceP
     return nothing
 end
 
-function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::SpaceParm{N,M,D}) where {T,N,M,D}
+function krnl_force_impr_pln!(frc1, frc2, U::AbstractArray{T}, c0, c1, ipl, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {T,N,M,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
diff --git a/src/YM/YMfields.jl b/src/YM/YMfields.jl
index dfcbbf2..df8094a 100644
--- a/src/YM/YMfields.jl
+++ b/src/YM/YMfields.jl
@@ -34,7 +34,7 @@ function randomize!(f, lp::SpaceParm, ymws::YMworkspace)
     return nothing
 end
 
-function krnl_assign_SU3!(frc, m, lp::SpaceParm)
+function krnl_assign_SU3!(frc, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {N,M,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
     for id in 1:lp.ndim
@@ -45,7 +45,7 @@ function krnl_assign_SU3!(frc, m, lp::SpaceParm)
     return nothing
 end
 
-function krnl_assign_SU2!(frc, m, lp::SpaceParm)
+function krnl_assign_SU2!(frc, m, lp::SpaceParm{N,M,BC_PERIODIC,D}) where {N,M,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
     for id in 1:lp.ndim
diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl
index c37fd37..5d4f959 100644
--- a/src/YM/YMflow.jl
+++ b/src/YM/YMflow.jl
@@ -25,7 +25,7 @@ function add_zth_term(ymws::YMworkspace, U, lp)
     return nothing
 end
 
-function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::SpaceParm{N,M,D}) where {TA,TG,N,M,D}
+function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::SpaceParm{N,M,B,D}) where {TA,TG,N,M,B,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -129,7 +129,7 @@ zfl_rk3(U, ns, eps, lp::SpaceParm, ymws::YMworkspace) = flw_rk3(U, ns, eps, 5.0/
 Measure the action density `E(t)` using the plaquette discretization. If the argument `Eslc`
 the contribution for each Euclidean time slice and plane are returned.
 """
-function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,D}, ymws::YMworkspace) where {T,N,M,D}
+function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) where {T,N,M,B,D}
 
     @timeit "E(t) plaquette measurement" begin
 
@@ -163,10 +163,10 @@ function Eoft_plaq(Eslc, U, gp::GaugeParm{T}, lp::SpaceParm{N,M,D}, ymws::YMwork
     return sum(Eslc)/lp.iL[end]
 end
 
-Eoft_plaq(U, gp::GaugeParm{T}, lp::SpaceParm{N,M,D}, ymws::YMworkspace) where {T,N,M,D} = Eoft_plaq(zeros(T,lp.iL[end],M), U, gp, lp, ymws)
+Eoft_plaq(U, gp::GaugeParm{T}, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace) where {T,N,M,B,D} = Eoft_plaq(zeros(T,lp.iL[end],M), U, gp, lp, ymws)
 
 
-function krnl_plaq_pln!(plx, U::AbstractArray{T}, ipl, lp::SpaceParm{N,M,D}) where {T,N,M,D}
+function krnl_plaq_pln!(plx, U::AbstractArray{T}, ipl, lp::SpaceParm{N,M,B,D}) where {T,N,M,B,D}
     
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -187,7 +187,7 @@ end
 Measure the topological charge `Q` of the configuration `U`. If the argument `Qslc` is present
 the contribution for each Euclidean time slice are returned.
 """
-function Qtop(Qslc, U, lp::SpaceParm{4,M,D}, ymws::YMworkspace) where {M,D}
+function Qtop(Qslc, U, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace) where {M,B,D}
 
     @timeit "Qtop measurement" begin
 
@@ -228,7 +228,7 @@ Qtop(U, lp::SpaceParm{4,M,D}, ymws::YMworkspace{T}) where {T,M,D} = Qtop(zeros(T
 Measure the action density `E(t)` using the clover discretization. If the argument `Eslc`
 the contribution for each Euclidean time slice and plane are returned.
 """
-function Eoft_clover(Eslc, U, lp::SpaceParm{4,M,D}, ymws::YMworkspace{T}) where {T,M,D}
+function Eoft_clover(Eslc, U, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace{T}) where {T,M,B,D}
 
     function acum(ipl1, ipl2, Etmp)
 
@@ -280,9 +280,9 @@ function Eoft_clover(Eslc, U, lp::SpaceParm{4,M,D}, ymws::YMworkspace{T}) where
 
     return sum(Eslc)/lp.iL[end]
 end
-Eoft_clover(U, lp::SpaceParm{N,M,D}, ymws::YMworkspace{T}) where {T,N,M,D} = Eoft_clover(zeros(T,lp.iL[end],M), U, lp, ymws)
+Eoft_clover(U, lp::SpaceParm{N,M,B,D}, ymws::YMworkspace{T}) where {T,N,M,B,D} = Eoft_clover(zeros(T,lp.iL[end],M), U, lp, ymws)
 
-function krnl_add_et!(rm, op, frc1, U, lp::SpaceParm{4,M,D}) where {M,D}
+function krnl_add_et!(rm, op, frc1, U, lp::SpaceParm{4,M,B,D}) where {M,B,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -294,7 +294,7 @@ function krnl_add_et!(rm, op, frc1, U, lp::SpaceParm{4,M,D}) where {M,D}
     return nothing
 end
 
-function krnl_add_qd!(rm, op, frc1, frc2, U, lp::SpaceParm{4,M,D}) where {M,D}
+function krnl_add_qd!(rm, op, frc1, frc2, U, lp::SpaceParm{4,M,B,D}) where {M,B,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
 
@@ -304,7 +304,7 @@ function krnl_add_qd!(rm, op, frc1, frc2, U, lp::SpaceParm{4,M,D}) where {M,D}
     return nothing
 end
 
-function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::SpaceParm{4,M,D}) where {T,M,D}
+function krnl_field_tensor!(frc1, frc2, U::AbstractArray{T}, ipl1, ipl2, lp::SpaceParm{4,M,B,D}) where {T,M,B,D}
 
     b, r = CUDA.threadIdx().x, CUDA.blockIdx().x