diff --git a/Project.toml b/Project.toml
index bbe490e4c..474fb968a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -34,7 +34,7 @@ LoggingExtras = "1"
 MPSKit = "0.12"
 MPSKitModels = "0.4"
 OhMyThreads = "0.7"
-OptimKit = "0.3"
+OptimKit = "0.3, 0.4"
 Printf = "1"
 QuadGK = "2.11.1"
 Random = "1"
diff --git a/README.md b/README.md
index 25f9a9248..381cb8d90 100644
--- a/README.md
+++ b/README.md
@@ -47,7 +47,7 @@ chi = 20
 ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi))
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
-    optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),
+    optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=3),
     gradient_alg=LinSolver(),
     reuse_env=true,
 )
diff --git a/docs/src/index.md b/docs/src/index.md
index 709f42222..115349e12 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -29,7 +29,7 @@ chi = 20
 ctm_alg = SimultaneousCTMRG(; tol=1e-10, trscheme=truncdim(chi))
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
-    optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),
+    optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=3),
     gradient_alg=LinSolver(),
     reuse_env=true,
 )
diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl
index 34baed189..3b0b0d657 100644
--- a/examples/heisenberg.jl
+++ b/examples/heisenberg.jl
@@ -14,7 +14,7 @@ H = heisenberg_XYZ(InfiniteSquare(); Jx=-1, Jy=1, Jz=-1)
 ctm_alg = SimultaneousCTMRG(; tol=1e-10, verbosity=2)
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
-    optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=2),
+    optimizer=LBFGS(4; maxiter=100, gradtol=1e-4, verbosity=3),
     gradient_alg=LinSolver(; solver=GMRES(; tol=1e-6, maxiter=100)),
     reuse_env=true,
 )
diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index b9ea58e64..314c80833 100644
--- a/src/PEPSKit.jl
+++ b/src/PEPSKit.jl
@@ -76,7 +76,7 @@ include("utility/symmetrization.jl")
         const ctmrg_alg = SimultaneousCTMRG(
             ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
         )
-        const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
+        const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=3)
         const gradient_linsolver = KrylovKit.BiCGStab(;
             maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
         )
@@ -130,7 +130,7 @@ module Defaults
     const ctmrg_alg = SimultaneousCTMRG(
         ctmrg_tol, ctmrg_maxiter, ctmrg_miniter, 2, projector_alg
     )
-    const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
+    const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=3)
     const gradient_linsolver = KrylovKit.BiCGStab(;
         maxiter=Defaults.fpgrad_maxiter, tol=Defaults.fpgrad_tol
     )
diff --git a/test/heisenberg.jl b/test/heisenberg.jl
index 29a0f0d5e..48e58be01 100644
--- a/test/heisenberg.jl
+++ b/test/heisenberg.jl
@@ -11,7 +11,7 @@ Dbond = 2
 χenv = 16
 ctm_alg = SimultaneousCTMRG()
 opt_alg = PEPSOptimize(;
-    boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2)
+    boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=3)
 )
 # compare against Juraj Hasik's data:
 # https://github.com/jurajHasik/j1j2_ipeps_states/blob/main/single-site_pg-C4v-A1/j20.0/state_1s_A1_j20.0_D2_chi_opt48.dat
diff --git a/test/j1j2_model.jl b/test/j1j2_model.jl
index 8dd883972..5d13bba74 100644
--- a/test/j1j2_model.jl
+++ b/test/j1j2_model.jl
@@ -11,7 +11,7 @@ using OptimKit
 ctm_alg = SimultaneousCTMRG()
 opt_alg = PEPSOptimize(;
     boundary_alg=ctm_alg,
-    optimizer=LBFGS(4; gradtol=1e-3, verbosity=2),
+    optimizer=LBFGS(4; gradtol=1e-3, verbosity=3),
     gradient_alg=LinSolver(; iterscheme=:diffgauge),
 )
 
diff --git a/test/pwave.jl b/test/pwave.jl
index dcdecf7ea..4e36b14e9 100644
--- a/test/pwave.jl
+++ b/test/pwave.jl
@@ -12,7 +12,7 @@ H = pwave_superconductor(InfiniteSquare(unitcell...))
 χenv = 16
 ctm_alg = SimultaneousCTMRG(; maxiter=150)
 opt_alg = PEPSOptimize(;
-    boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=10, gradtol=1e-3, verbosity=2)
+    boundary_alg=ctm_alg, optimizer=LBFGS(4; maxiter=10, gradtol=1e-3, verbosity=3)
 )
 
 # initialize states
diff --git a/test/tf_ising.jl b/test/tf_ising.jl
index 8b84696f2..903770214 100644
--- a/test/tf_ising.jl
+++ b/test/tf_ising.jl
@@ -21,12 +21,12 @@ mˣ = 0.91
 χenv = 16
 ctm_alg = SimultaneousCTMRG()
 opt_alg = PEPSOptimize(;
-    boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=2)
+    boundary_alg=ctm_alg, optimizer=LBFGS(4; gradtol=1e-3, verbosity=3)
 )
 
 # initialize states
 H = transverse_field_ising(InfiniteSquare(); g)
-Random.seed!(91283219347)
+Random.seed!(2928528935)
 psi_init = InfinitePEPS(2, χbond)
 env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, ctm_alg)