diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index 941657a0..75588875 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -33,14 +33,29 @@ jobs:
         with:
           version: ${{ matrix.version }}
 
+      - name: Set up Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: '3.10'  # Specify your Python version
+
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip 
+          pip install torch
+          pip install numpy
+
       - name: Cache Julia artifacts
         uses: julia-actions/cache@v2
 
       - name: Build project
         uses: julia-actions/julia-buildpkg@v1
+        env:
+          PYTHON: python
 
       - name: Run tests
         uses: julia-actions/julia-runtest@v1
+        with:
+          test_args: "run-nn-opinf-tests"
 
       - name: Process coverage
         uses: julia-actions/julia-processcoverage@v1
diff --git a/Project.toml b/Project.toml
index 345650cd..02ed14d8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,6 +20,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
+PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 
 [compat]
 DelimitedFiles = "1"
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1-for-test.yaml b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1-for-test.yaml
new file mode 100644
index 00000000..ccc396e4
--- /dev/null
+++ b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1-for-test.yaml
@@ -0,0 +1,40 @@
+type: single
+input mesh file: cuboid-1.g
+output mesh file: cuboid-1.e
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      fine: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+  Schwarz overlap:
+    - side set: ssz+
+      source: cuboid-2-for-test.yaml
+      source block: coarse
+      source side set: ssz-
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1.yaml b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1.yaml
new file mode 100644
index 00000000..8e56c8f2
--- /dev/null
+++ b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1.yaml
@@ -0,0 +1,40 @@
+type: single
+input mesh file: cuboid-1.g
+output mesh file: cuboid-1.e
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      fine: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+  Schwarz overlap:
+    - side set: ssz+
+      source: cuboid-2.yaml
+      source block: coarse
+      source side set: ssz-
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2-for-test.yaml b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2-for-test.yaml
new file mode 100644
index 00000000..84ea9cb9
--- /dev/null
+++ b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2-for-test.yaml
@@ -0,0 +1,41 @@
+type: single
+input mesh file: cuboid-2.g
+output mesh file: cuboid-2.e
+model:
+  type: neural network opinf rom
+  model-directory: neural-network-model/
+  ensemble-size: 0
+  material:
+    blocks:
+      coarse: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  OpInf Dirichlet:
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+      model-directory: neural-network-model/ 
+      ensemble-size: 0
+
+  OpInf Schwarz overlap:
+    - side set: ssz- 
+      source: cuboid-1-for-test.yaml
+      source block: fine 
+      source side set: ssz+
+      model-directory: neural-network-model/ 
+      ensemble-size: 0
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 1
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2.yaml b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2.yaml
new file mode 100644
index 00000000..a3d84e7a
--- /dev/null
+++ b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2.yaml
@@ -0,0 +1,41 @@
+type: single
+input mesh file: cuboid-2.g
+output mesh file: cuboid-2.e
+model:
+  type: neural network opinf rom
+  model-directory: neural-network-model/
+  ensemble-size: 5
+  material:
+    blocks:
+      coarse: hyperelastic
+    hyperelastic:
+      model: neohookean
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: Newmark
+  β: 0.25
+  γ: 0.5
+boundary conditions:
+  OpInf Dirichlet:
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+      model-directory: neural-network-model/ 
+      ensemble-size: 5
+
+  OpInf Schwarz overlap:
+    - side set: ssz- 
+      source: cuboid-1.yaml
+      source block: fine 
+      source side set: ssz+
+      model-directory: neural-network-model/ 
+      ensemble-size: 5
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-10
+  absolute tolerance: 1.0e-06
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids-for-test.yaml b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids-for-test.yaml
new file mode 100644
index 00000000..c67bbe5b
--- /dev/null
+++ b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids-for-test.yaml
@@ -0,0 +1,11 @@
+type: multi
+domains: ["cuboid-1-for-test.yaml", "cuboid-2-for-test.yaml"]
+Exodus output interval: 1
+CSV output interval: 1
+initial time: 0.0
+final time: 0.01
+time step: 0.01
+minimum iterations: 1
+maximum iterations: 1
+relative tolerance: 1.0e-12
+absolute tolerance: 1.0e-08
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids.yaml b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids.yaml
new file mode 100644
index 00000000..f5fd5338
--- /dev/null
+++ b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids.yaml
@@ -0,0 +1,11 @@
+type: multi
+domains: ["cuboid-1.yaml", "cuboid-2.yaml"]
+Exodus output interval: 1
+CSV output interval: 1
+initial time: 0.0
+final time: 1.0
+time step: 0.01
+minimum iterations: 1
+maximum iterations: 16
+relative tolerance: 1.0e-12
+absolute tolerance: 1.0e-08
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-0.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-0.pt
new file mode 100644
index 00000000..27d3778f
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-0.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-1.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-1.pt
new file mode 100644
index 00000000..852b8c35
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-1.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-2.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-2.pt
new file mode 100644
index 00000000..f0383fbc
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-2.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-3.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-3.pt
new file mode 100644
index 00000000..87cfa59a
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-3.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-4.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-4.pt
new file mode 100644
index 00000000..8cdc7d6d
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-nsz+-z-4.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--0.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--0.pt
new file mode 100644
index 00000000..0e8ee841
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--0.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--1.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--1.pt
new file mode 100644
index 00000000..64acf53d
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--1.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--2.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--2.pt
new file mode 100644
index 00000000..861b8acd
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--2.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--3.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--3.pt
new file mode 100644
index 00000000..a233566d
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--3.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--4.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--4.pt
new file mode 100644
index 00000000..390bf7b8
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/BC-ssz--4.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis-nsz+-z.npz b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis-nsz+-z.npz
new file mode 100644
index 00000000..5d3901b3
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis-nsz+-z.npz differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis-ssz-.npz b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis-ssz-.npz
new file mode 100644
index 00000000..d679ad19
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis-ssz-.npz differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis.npz b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis.npz
new file mode 100644
index 00000000..dd005f4b
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/nn-opinf-basis.npz differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model.pt
new file mode 100644
index 00000000..68cf51d7
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model_not_scaled.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model_not_scaled.pt
new file mode 100644
index 00000000..e3ff6de4
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model_not_scaled.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model_training_stats.npz b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model_training_stats.npz
new file mode 100644
index 00000000..dace0767
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/pytorch-model_training_stats.npz differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-0.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-0.pt
new file mode 100644
index 00000000..73cf8693
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-0.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-1.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-1.pt
new file mode 100644
index 00000000..edc68d04
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-1.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-2.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-2.pt
new file mode 100644
index 00000000..308aeff5
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-2.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-3.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-3.pt
new file mode 100644
index 00000000..770e090d
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-3.pt differ
diff --git a/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-4.pt b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-4.pt
new file mode 100644
index 00000000..251e54a2
Binary files /dev/null and b/examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model/stiffness-4.pt differ
diff --git a/examples/ahead/overlap/torsion/dynamic-opinf-fom/make_op_inf_models.py b/examples/ahead/overlap/torsion/dynamic-opinf-fom/make_op_inf_models.py
index 8e154264..40ca45e5 100644
--- a/examples/ahead/overlap/torsion/dynamic-opinf-fom/make_op_inf_models.py
+++ b/examples/ahead/overlap/torsion/dynamic-opinf-fom/make_op_inf_models.py
@@ -1,117 +1,23 @@
-import numpy as np
-import opinf
-import os
-from matplotlib import pyplot as plt
 import normaopinf
-import normaopinf.readers
-import normaopinf.calculus
-import romtools
-
-if __name__ == "__main__":
-    # Load in snapshots
-    cur_dir = os.getcwd()
-    solution_id = 2
-    displacement_snapshots,times = normaopinf.readers.load_displacement_csv_files(solution_directory=cur_dir,solution_id=solution_id,skip_files=1)
-    
-    # Identify which DOFs are free
-    free_dofs = normaopinf.readers.get_free_dofs(solution_directory=cur_dir,solution_id=solution_id)
-    
-    # Set values = 0 if DOFs are fixed
-    displacement_snapshots[free_dofs[:,:]==False] = 0.
-    
-    #Get sideset snapshots
-    sidesets = ["-Z_ss"]
-    sideset_snapshots = normaopinf.readers.load_sideset_displacement_csv_files(solution_directory=cur_dir,sidesets=sidesets,solution_id=solution_id,skip_files=1)
-    
-    # Create an energy-based truncater 
-    #tolerance = 1.e-5
-    #my_energy_truncater = romtools.vector_space.utils.EnergyBasedTruncater(1. - tolerance)
-    my_energy_truncater = romtools.vector_space.utils.BasisSizeTruncater(4)
-
-    # Now load in sidesets and create reduced spaces
-    # Note that I construct a separate basis for each x,y,z component. This isn't necessary 
-    ss_tspace = {}
-    reduced_sideset_snapshots = {}
-    for sideset in sidesets:
-        if sideset_snapshots[sideset].shape[0] ==  1:
-          ss_tspace[sideset] = romtools.VectorSpaceFromPOD(snapshots=sideset_snapshots[sideset],
-                                              truncater=my_energy_truncater,
-                                              shifter = None,
-                                              orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
-                                              scaler = romtools.vector_space.utils.NoOpScaler())
-          # Compute L2 orthogonal projection onto trial spaces
-          reduced_sideset_snapshots[sideset] = romtools.rom.optimal_l2_projection(sideset_snapshots[sideset],ss_tspace[sideset]) 
-        else:
-          comp_trial_space = []
-          for i in range(0,3):
-            tspace = romtools.VectorSpaceFromPOD(snapshots=sideset_snapshots[sideset][i:i+1],
-                                              truncater=my_energy_truncater,
-                                              shifter = None,
-                                              orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
-                                              scaler = romtools.vector_space.utils.NoOpScaler())
-            comp_trial_space.append(tspace)
-          ss_tspace[sideset] = romtools.CompositeVectorSpace(comp_trial_space)
-          reduced_sideset_snapshots[sideset] = romtools.rom.optimal_l2_projection(sideset_snapshots[sideset],ss_tspace[sideset]) 
-
-
-    ## Stack sidesets into one matrix for OpInf 
-    reduced_stacked_sideset_snapshots = None
-    for sideset in sidesets:
-        if reduced_stacked_sideset_snapshots is None:
-            reduced_stacked_sideset_snapshots = reduced_sideset_snapshots[sideset]*1.
-        else: 
-            reduced_stacked_sideset_snapshots = np.append(reduced_stacked_sideset_snapshots,reduced_sideset_snapshots[sideset],axis=0)
-    
-    # Create trial space for displacement vector
-    # Note again that I construct a separate basis for each x,y,z component. This isn't necessary 
-    trial_spaces = []
-    
-    my_energy_truncater = romtools.vector_space.utils.BasisSizeTruncater(60)
-
-    for i in range(0,3):
-      trial_space = romtools.VectorSpaceFromPOD(snapshots=displacement_snapshots[i:i+1],
-                                              truncater=my_energy_truncater,
-                                              shifter = None,
-                                              orthogonalizer=romtools.vector_space.utils.EuclideanL2Orthogonalizer(),
-                                              scaler = romtools.vector_space.utils.NoOpScaler())
-      trial_spaces.append(trial_space)
-    
-    trial_space = romtools.CompositeVectorSpace(trial_spaces)
-    
-    # Compute L2 orthogonal projection onto trial spaces 
-    uhat = romtools.rom.optimal_l2_projection(displacement_snapshots,trial_space)
-    u_ddots,times_dummy = normaopinf.readers.load_acceleration_csv_files(solution_directory=cur_dir,solution_id=solution_id,skip_files=1)
-    # Set values = 0 if DOFs are fixed
-    u_ddots[free_dofs[:,:]==False] = 0.
-    uhat_ddots = romtools.rom.optimal_l2_projection(u_ddots*1.,trial_space)
-   
-    # Construct an opinf "AB" model (linear in the state and linear in the exogenous inputs)
-    #   Note: I don't construct a cAB ROM in this example since I know there is no forcing vector
-    l2solver = opinf.lstsq.L2Solver(regularizer=1e-11)
-    opinf_model = opinf.models.ContinuousModel("AB",solver=l2solver)
-    opinf_model.fit(states=uhat, ddts=uhat_ddots,inputs=reduced_stacked_sideset_snapshots)
-
-    # Flip signs to match convention of K on the LHS
-    K = -opinf_model.A_.entries
-    B = opinf_model.B_.entries
-
-    ## Now extract boundary operators and create dictionary to save
-    col_start = 0
-    sideset_operators = {}
-    for sideset in sidesets:
-      num_dofs = reduced_sideset_snapshots[sideset].shape[0]
-      val = np.einsum('kr,vnr->vkn',B[:,col_start:col_start + num_dofs] , ss_tspace[sideset].get_basis() )
-      shape2 = B[:,col_start:col_start + num_dofs] @ ss_tspace[sideset].get_basis()[0].transpose()
-      sideset_operators["B_" + sideset] = val#
-      col_start += num_dofs
-    
-    f = np.zeros(K.shape[0])
-    vals_to_save = sideset_operators 
-    vals_to_save["basis"] = trial_space.get_basis() 
-    vals_to_save["K"] = K
-    vals_to_save["f"] = f
-    
-    np.savez('opinf-operator',**vals_to_save)
-
-
+import normaopinf.opinf
+import os
+import numpy as np
 
+if __name__ == '__main__':
+  settings = {}
+  settings['fom-yaml-file'] = "torsion.yaml"
+  settings['training-data-directories'] = [os.getcwd()]
+  settings['solution-id'] = 1
+  settings['model-type'] = 'quadratic'
+  settings['forcing'] =  False
+  settings['truncation-type'] = 'size'
+  settings['boundary-truncation-type'] =  'energy' 
+  settings['regularization-parameter'] =  'automatic'
+  settings['trial-space-splitting-type'] = 'split'
+  settings['acceleration-computation-type'] = 'finite-difference'
+  sizes = np.array([2,4,6,8,10,12,14,16,18,20,40,60],dtype=int)
+  for i,size in enumerate(sizes):
+    settings['truncation-value'] = int(size)
+    settings['boundary-truncation-value'] = 1. - 1.e-5 
+    settings['operator-name'] = 'quadratic-opinf-operator-rom-dim-' + str(size)
+    normaopinf.opinf.make_opinf_model(settings)
diff --git a/src/ics_bcs.jl b/src/ics_bcs.jl
index 2d4ba835..1c422ea1 100644
--- a/src/ics_bcs.jl
+++ b/src/ics_bcs.jl
@@ -3,6 +3,8 @@
 # the U.S. Government retains certain rights in this software. This software
 # is released under the BSD license detailed in the file license.txt in the
 # top-level Norma.jl directory.
+using PyCall
+import NPZ
 
 @variables t, x, y, z
 D = Differential(t)
@@ -242,32 +244,6 @@ function SMCouplingSchwarzBC(
     end
 end
 
-function apply_bc(model::OpInfModel, bc::SMDirichletBC)
-    model.fom_model.time = model.time
-    apply_bc(model.fom_model, bc)
-    bc_vector = zeros(0)
-    for node_index in bc.node_set_node_indices
-        disp_val = model.fom_model.current[bc.offset, node_index] - model.fom_model.reference[bc.offset, node_index]
-        push!(bc_vector, disp_val)
-    end
-
-    offset = bc.offset
-    if offset == 1
-        offset_name = "x"
-    end
-    if offset == 2
-        offset_name = "y"
-    end
-    if offset == 3
-        offset_name = "z"
-    end
-
-    op_name = "B_" * bc.node_set_name * "-" * offset_name
-    bc_operator = model.opinf_rom[op_name]
-    # SM Dirichlet BC are only defined on a single x,y,z
-    return model.reduced_boundary_forcing[:] += bc_operator[1, :, :] * bc_vector
-end
-
 function apply_bc(model::SolidMechanics, bc::SMDirichletBC)
     for node_index in bc.node_set_node_indices
         values = Dict(
@@ -404,10 +380,6 @@ function find_point_in_mesh(point::Vector{Float64}, model::SolidMechanics, blk_i
     return node_indices, ξ, found
 end
 
-function find_point_in_mesh(point::Vector{Float64}, model::RomModel, blk_id::Int, tol::Float64)
-    node_indices, ξ, found = find_point_in_mesh(point, model.fom_model, blk_id, tol)
-    return node_indices, ξ, found
-end
 
 function apply_bc_detail(model::SolidMechanics, bc::SMContactSchwarzBC)
     if bc.is_dirichlet == true
@@ -425,27 +397,6 @@ function apply_bc_detail(model::SolidMechanics, bc::CouplingSchwarzBoundaryCondi
     end
 end
 
-function apply_bc_detail(model::OpInfModel, bc::CouplingSchwarzBoundaryCondition)
-    if bc.coupled_subsim.model isa SolidMechanics
-        ## Apply BC to the FOM vector
-        apply_bc_detail(model.fom_model, bc)
-
-        # populate our own BC vector
-        bc_vector = zeros(3, length(bc.side_set_node_indices))
-        for i in 1:length(bc.side_set_node_indices)
-            node_index = bc.side_set_node_indices[i]
-            bc_vector[:, i] = model.fom_model.current[:, node_index] - model.fom_model.reference[:, node_index]
-        end
-        op_name = "B_" * bc.side_set_name
-        bc_operator = model.opinf_rom[op_name]
-        for i in 1:3
-            model.reduced_boundary_forcing[:] += bc_operator[i, :, :] * bc_vector[i, :]
-        end
-    else
-        throw("ROM-ROM coupling not supported yet")
-    end
-end
-
 function apply_sm_schwarz_coupling_dirichlet(model::SolidMechanics, bc::CouplingSchwarzBoundaryCondition)
     if bc.coupled_subsim.model isa SolidMechanics
         for i in 1:length(bc.side_set_node_indices)
@@ -560,78 +511,6 @@ function apply_bc(model::SolidMechanics, bc::SchwarzBoundaryCondition)
     return copy_solution_source_targets(bc.coupled_subsim.integrator, bc.coupled_subsim.solver, bc.coupled_subsim.model)
 end
 
-function apply_bc(model::RomModel, bc::SchwarzBoundaryCondition)
-    global_sim = bc.coupled_subsim.params["global_simulation"]
-    controller = global_sim.controller
-    if bc isa SMContactSchwarzBC && controller.active_contact == false
-        return nothing
-    end
-    empty_history = length(controller.time_hist) == 0
-    same_step = controller.same_step
-    if empty_history == true
-        apply_bc_detail(model, bc)
-        return nothing
-    end
-    # Save solution of coupled simulation
-    saved_disp = bc.coupled_subsim.integrator.displacement
-    saved_velo = bc.coupled_subsim.integrator.velocity
-    saved_acce = bc.coupled_subsim.integrator.acceleration
-    saved_∂Ω_f = bc.coupled_subsim.model.internal_force
-    time = model.time
-    coupled_name = bc.coupled_subsim.name
-    coupled_index = global_sim.subsim_name_index_map[coupled_name]
-    time_hist = controller.time_hist[coupled_index]
-    disp_hist = controller.disp_hist[coupled_index]
-    velo_hist = controller.velo_hist[coupled_index]
-    acce_hist = controller.acce_hist[coupled_index]
-    ∂Ω_f_hist = controller.∂Ω_f_hist[coupled_index]
-    interp_disp = same_step == true ? disp_hist[end] : interpolate(time_hist, disp_hist, time)
-    interp_velo = same_step == true ? velo_hist[end] : interpolate(time_hist, velo_hist, time)
-    interp_acce = same_step == true ? acce_hist[end] : interpolate(time_hist, acce_hist, time)
-    interp_∂Ω_f = same_step == true ? ∂Ω_f_hist[end] : interpolate(time_hist, ∂Ω_f_hist, time)
-    if bc.coupled_subsim.model isa SolidMechanics
-        bc.coupled_subsim.model.internal_force = interp_∂Ω_f
-    elseif bc.coupled_subsim.model isa RomModel
-        bc.coupled_subsim.model.fom_model.internal_force = interp_∂Ω_f
-    end
-
-    if bc isa SMContactSchwarzBC || bc isa SMNonOverlapSchwarzBC
-        relaxation_parameter = controller.relaxation_parameter
-        schwarz_iteration = controller.iteration_number
-        if schwarz_iteration == 1
-            lambda_dispᵖʳᵉᵛ = zeros(length(interp_disp))
-            lambda_veloᵖʳᵉᵛ = zeros(length(interp_velo))
-            lambda_acceᵖʳᵉᵛ = zeros(length(interp_acce))
-        else
-            lambda_dispᵖʳᵉᵛ = controller.lambda_disp[coupled_index]
-            lambda_veloᵖʳᵉᵛ = controller.lambda_velo[coupled_index]
-            lambda_acceᵖʳᵉᵛ = controller.lambda_acce[coupled_index]
-        end
-        bc.coupled_subsim.integrator.displacement =
-            controller.lambda_disp[coupled_index] =
-                relaxation_parameter * interp_disp + (1 - relaxation_parameter) * lambda_dispᵖʳᵉᵛ
-        bc.coupled_subsim.integrator.velocity =
-            controller.lambda_velo[coupled_index] =
-                relaxation_parameter * interp_velo + (1 - relaxation_parameter) * lambda_veloᵖʳᵉᵛ
-        bc.coupled_subsim.integrator.acceleration =
-            controller.lambda_acce[coupled_index] =
-                relaxation_parameter * interp_acce + (1 - relaxation_parameter) * lambda_acceᵖʳᵉᵛ
-    else
-        bc.coupled_subsim.integrator.displacement = interp_disp
-        bc.coupled_subsim.integrator.velocity = interp_velo
-        bc.coupled_subsim.integrator.acceleration = interp_acce
-    end
-    # Copies from integrator to model
-    copy_solution_source_targets(bc.coupled_subsim.integrator, bc.coupled_subsim.solver, bc.coupled_subsim.model)
-    apply_bc_detail(model, bc)
-    bc.coupled_subsim.integrator.displacement = saved_disp
-    bc.coupled_subsim.integrator.velocity = saved_velo
-    bc.coupled_subsim.integrator.acceleration = saved_acce
-    bc.coupled_subsim.model.internal_force = saved_∂Ω_f
-    # Copy from integrator to model
-    return copy_solution_source_targets(bc.coupled_subsim.integrator, bc.coupled_subsim.solver, bc.coupled_subsim.model)
-end
-
 function transfer_normal_component(source::Vector{Float64}, target::Vector{Float64}, normal::Vector{Float64})
     normal_projection = normal * normal'
     tangent_projection = I(length(normal)) - normal_projection
@@ -834,6 +713,9 @@ function create_bcs(params::Parameters)
             if bc_type == "Dirichlet"
                 boundary_condition = SMDirichletBC(input_mesh, bc_setting_params)
                 push!(boundary_conditions, boundary_condition)
+            elseif bc_type == "OpInf Dirichlet" 
+                boundary_condition = SMOpInfDirichletBC(input_mesh, bc_setting_params)
+                push!(boundary_conditions, boundary_condition)
             elseif bc_type == "Neumann"
                 boundary_condition = SMNeumannBC(input_mesh, bc_setting_params)
                 push!(boundary_conditions, boundary_condition)
@@ -859,6 +741,17 @@ function create_bcs(params::Parameters)
                 coupled_subsim = sim.subsims[coupled_subdomain_index]
                 boundary_condition = SMCouplingSchwarzBC(subsim, coupled_subsim, input_mesh, bc_type, bc_setting_params)
                 push!(boundary_conditions, boundary_condition)
+            elseif bc_type == "OpInf Schwarz overlap"
+                sim = params["global_simulation"]
+                subsim_name = params["name"]
+                subdomain_index = sim.subsim_name_index_map[subsim_name]
+                subsim = sim.subsims[subdomain_index]
+                coupled_subsim_name = bc_setting_params["source"]
+                coupled_subdomain_index = sim.subsim_name_index_map[coupled_subsim_name]
+                coupled_subsim = sim.subsims[coupled_subdomain_index]
+                boundary_condition = SMOpInfCouplingSchwarzBC(subsim, coupled_subsim, input_mesh, bc_type, bc_setting_params)
+                push!(boundary_conditions, boundary_condition)
+
             else
                 error("Unknown boundary condition type : ", bc_type)
             end
@@ -880,13 +773,6 @@ function apply_bcs(model::SolidMechanics)
     end
 end
 
-function apply_bcs(model::RomModel)
-    model.reduced_boundary_forcing[:] .= 0.0
-    for boundary_condition in model.boundary_conditions
-        apply_bc(model, boundary_condition)
-    end
-end
-
 function assign_velocity!(
     velocity::Matrix{Float64}, offset::Int64, node_index::Int32, velo_val::Float64, context::String
 )
@@ -950,34 +836,6 @@ function apply_ics(params::Parameters, model::SolidMechanics)
     end
 end
 
-function apply_ics(params::Parameters, model::RomModel)
-    apply_ics(params, model.fom_model)
-
-    if haskey(params, "initial conditions") == false
-        return nothing
-    end
-    n_var, n_node, n_mode = model.basis.size
-    n_var_fom, n_node_fom = size(model.fom_model.current)
-
-    # Make sure basis is the right size
-    if n_var != n_var_fom || n_node != n_node_fom
-        throw("Basis is wrong size")
-    end
-
-    # project onto basis
-    for k in 1:n_mode
-        model.reduced_state[k] = 0.0
-        model.reduced_velocity[k] = 0.0
-        for j in 1:n_node
-            for n in 1:n_var
-                model.reduced_state[k] +=
-                    model.basis[n, j, k] * (model.fom_model.current[n, j] - model.fom_model.reference[n, j])
-                model.reduced_velocity[k] += model.basis[n, j, k] * (model.fom_model.velocity[n, j])
-            end
-        end
-    end
-end
-
 function pair_schwarz_bcs(sim::MultiDomainSimulation)
     for subsim in sim.subsims
         model = subsim.model
diff --git a/src/ics_bcs_types.jl b/src/ics_bcs_types.jl
index d450e14e..586a1d6e 100644
--- a/src/ics_bcs_types.jl
+++ b/src/ics_bcs_types.jl
@@ -26,6 +26,19 @@ mutable struct SMDirichletBC <: RegularBoundaryCondition
     acce_num::Num
 end
 
+mutable struct SMOpInfDirichletBC <: RegularBoundaryCondition 
+    node_set_name::String
+    offset::Int64
+    node_set_id::Int64
+    node_set_node_indices::Vector{Int64}
+    disp_num::Num
+    velo_num::Num
+    acce_num::Num
+    fom_bc::SMDirichletBC
+    nn_model::Any
+    basis::Any
+end
+
 mutable struct SMDirichletInclined <: RegularBoundaryCondition
     node_set_name::String
     node_set_id::Int64
@@ -74,6 +87,21 @@ mutable struct SMOverlapSchwarzBC <: OverlapSchwarzBoundaryCondition
     swap_bcs::Bool
 end
 
+mutable struct SMOpInfOverlapSchwarzBC <: OverlapSchwarzBoundaryCondition
+    side_set_name::String
+    side_set_node_indices::Vector{Int64}
+    coupled_nodes_indices::Vector{Vector{Int64}}
+    interpolation_function_values::Vector{Vector{Float64}}
+    coupled_subsim::Simulation
+    subsim::Simulation
+    is_dirichlet::Bool
+    swap_bcs::Bool
+    fom_bc::SMOverlapSchwarzBC
+    nn_model::Any
+    basis::Any
+end
+
+
 mutable struct SMNonOverlapSchwarzBC <: NonOverlapSchwarzBoundaryCondition
     side_set_id::Int64
     side_set_node_indices::Vector{Int64}
diff --git a/src/model.jl b/src/model.jl
index d0b23c08..56cbcb50 100644
--- a/src/model.jl
+++ b/src/model.jl
@@ -7,79 +7,8 @@
 include("constitutive.jl")
 include("interpolation.jl")
 include("ics_bcs.jl")
-
+include("opinf/opinf_ics_bcs.jl")
 using Base.Threads: @threads, threadid, nthreads
-using NPZ
-
-function LinearOpInfRom(params::Parameters)
-    params["mesh smoothing"] = false
-    fom_model = SolidMechanics(params)
-    reference = fom_model.reference
-    opinf_model_file = params["model"]["model-file"]
-    opinf_model = NPZ.npzread(opinf_model_file)
-    basis = opinf_model["basis"]
-    _, _, reduced_dim = size(basis)
-    num_dofs = reduced_dim
-    time = 0.0
-    failed = false
-    null_vec = zeros(num_dofs)
-
-    reduced_state = zeros(num_dofs)
-    reduced_velocity = zeros(num_dofs)
-    reduced_boundary_forcing = zeros(num_dofs)
-    free_dofs = trues(num_dofs)
-    boundary_conditions = Vector{BoundaryCondition}()
-    return LinearOpInfRom(
-        opinf_model,
-        basis,
-        reduced_state,
-        reduced_velocity,
-        reduced_boundary_forcing,
-        null_vec,
-        free_dofs,
-        boundary_conditions,
-        time,
-        failed,
-        fom_model,
-        reference,
-        false,
-    )
-end
-
-function QuadraticOpInfRom(params::Parameters)
-    params["mesh smoothing"] = false
-    fom_model = SolidMechanics(params)
-    reference = fom_model.reference
-    opinf_model_file = params["model"]["model-file"]
-    opinf_model = NPZ.npzread(opinf_model_file)
-    basis = opinf_model["basis"]
-    _, _, reduced_dim = size(basis)
-    num_dofs = reduced_dim
-    time = 0.0
-    failed = false
-    null_vec = zeros(num_dofs)
-
-    reduced_state = zeros(num_dofs)
-    reduced_velocity = zeros(num_dofs)
-    reduced_boundary_forcing = zeros(num_dofs)
-    free_dofs = trues(num_dofs)
-    boundary_conditions = Vector{BoundaryCondition}()
-    return QuadraticOpInfRom(
-        opinf_model,
-        basis,
-        reduced_state,
-        reduced_velocity,
-        reduced_boundary_forcing,
-        null_vec,
-        free_dofs,
-        boundary_conditions,
-        time,
-        failed,
-        fom_model,
-        reference,
-        false,
-    )
-end
 
 function SolidMechanics(params::Parameters)
     input_mesh = params["input_mesh"]
@@ -208,6 +137,8 @@ function create_model(params::Parameters)
         return LinearOpInfRom(params)
     elseif model_name == "quadratic opinf rom"
         return QuadraticOpInfRom(params)
+    elseif model_name == "neural network opinf rom"
+        return NeuralNetworkOpInfRom(params)
 
     else
         error("Unknown type of model : ", model_name)
@@ -779,4 +710,4 @@ function get_block_connectivity(mesh::ExodusDatabase, blk_id::Integer)
     _, num_elems, num_nodes, _, _, _ = Exodus.read_block_parameters(mesh, Int32(blk_id))
     conn = Exodus.read_block_connectivity(mesh, Int32(blk_id), num_elems * num_nodes)
     return reshape(conn, (num_elems, num_nodes))
-end
\ No newline at end of file
+end
diff --git a/src/model_types.jl b/src/model_types.jl
index 391ed8fe..6cbf6919 100644
--- a/src/model_types.jl
+++ b/src/model_types.jl
@@ -5,8 +5,6 @@
 # top-level Norma.jl directory.
 
 abstract type Model end
-abstract type RomModel <: Model end
-abstract type OpInfModel <: RomModel end
 using SparseArrays
 
 @enum Kinematics begin
@@ -95,36 +93,3 @@ mutable struct SolidMechanics <: Model
     kinematics::Kinematics
 end
 
-mutable struct QuadraticOpInfRom <: OpInfModel
-    opinf_rom::Dict{Any,Any}
-    basis::Array{Float64}
-    reduced_state::Vector{Float64}
-    reduced_velocity::Vector{Float64}
-    reduced_boundary_forcing::Vector{Float64}
-    #internal_force not used, but include to ease interfacing in Schwarz
-    internal_force::Vector{Float64}
-    free_dofs::BitVector
-    boundary_conditions::Vector{BoundaryCondition}
-    time::Float64
-    failed::Bool
-    fom_model::SolidMechanics
-    reference::Matrix{Float64}
-    inclined_support::Bool
-end
-
-mutable struct LinearOpInfRom <: OpInfModel
-    opinf_rom::Dict{Any,Any}
-    basis::Array{Float64}
-    reduced_state::Vector{Float64}
-    reduced_velocity::Vector{Float64}
-    reduced_boundary_forcing::Vector{Float64}
-    #internal_force not used, but include to ease interfacing in Schwarz
-    internal_force::Vector{Float64}
-    free_dofs::BitVector
-    boundary_conditions::Vector{BoundaryCondition}
-    time::Float64
-    failed::Bool
-    fom_model::SolidMechanics
-    reference::Matrix{Float64}
-    inclined_support::Bool
-end
diff --git a/src/opinf/opinf_ics_bcs.jl b/src/opinf/opinf_ics_bcs.jl
new file mode 100644
index 00000000..8fe45c0d
--- /dev/null
+++ b/src/opinf/opinf_ics_bcs.jl
@@ -0,0 +1,387 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+using PyCall
+import NPZ
+
+@variables t, x, y, z
+D = Differential(t)
+
+
+function SMOpInfDirichletBC(input_mesh::ExodusDatabase, bc_params::Dict{String,Any})
+    fom_bc = SMDirichletBC(input_mesh,bc_params)
+    node_set_name = bc_params["node set"]
+    expression = bc_params["function"]
+    offset = component_offset_from_string(bc_params["component"])
+    node_set_id = node_set_id_from_name(node_set_name, input_mesh)
+    node_set_node_indices = Exodus.read_node_set_nodes(input_mesh, node_set_id)
+    # expression is an arbitrary function of t, x, y, z in the input file
+    disp_num = eval(Meta.parse(expression))
+    velo_num = expand_derivatives(D(disp_num))
+    acce_num = expand_derivatives(D(velo_num))
+
+    opinf_model_directory = bc_params["model-directory"]
+    py""" 
+    import torch
+    def get_model(model_file):
+      import os
+      assert os.path.isfile(model_file) ,  print(model_file + " cannot be found" )
+      return torch.load(model_file)
+    """
+    ensemble_size = bc_params["ensemble-size"]
+
+    if offset == 1
+        offset_name = "x"
+    end
+    if offset == 2
+        offset_name = "y"
+    end
+    if offset == 3
+        offset_name = "z"
+    end
+
+
+    model = []
+    for i in 1:ensemble_size
+      tmp =  py"get_model"(opinf_model_directory * "/BC-" * node_set_name * "-" * offset_name * "-" * string(i-1) * ".pt")
+      push!(model,tmp)
+    end
+
+    basis_file = bc_params["model-directory"] * "/nn-opinf-basis-" * node_set_name * "-" * offset_name * ".npz" 
+    basis = NPZ.npzread(basis_file)
+    basis = basis["basis"]
+
+
+    SMOpInfDirichletBC(
+        node_set_name,
+        offset,
+        node_set_id,
+        node_set_node_indices,
+        disp_num,
+        velo_num,
+        acce_num,
+        fom_bc,
+        model,
+        basis,
+    )
+end
+
+
+function SMOpInfCouplingSchwarzBC(
+    subsim::SingleDomainSimulation,
+    coupled_subsim::SingleDomainSimulation,
+    input_mesh::ExodusDatabase,
+    bc_type::String,
+    bc_params::Parameters,
+)
+    fom_bc = SMCouplingSchwarzBC(subsim,coupled_subsim,input_mesh,"Schwarz overlap",bc_params) 
+    side_set_name = bc_params["side set"]
+    side_set_id = side_set_id_from_name(side_set_name, input_mesh)
+    _, _, side_set_node_indices = get_side_set_local_from_global_map(input_mesh, side_set_id)
+    coupled_block_name = bc_params["source block"]
+    if typeof(coupled_subsim.model) <: RomModel
+        coupled_mesh = coupled_subsim.model.fom_model.mesh
+    else
+        coupled_mesh = coupled_subsim.model.mesh
+    end
+    coupled_block_id = block_id_from_name(coupled_block_name, coupled_mesh)
+    element_type_string = Exodus.read_block_parameters(coupled_mesh, coupled_block_id)[1]
+    element_type = element_type_from_string(element_type_string)
+    coupled_side_set_name = bc_params["source side set"]
+    coupled_side_set_id = side_set_id_from_name(coupled_side_set_name, coupled_mesh)
+    coupled_nodes_indices = Vector{Vector{Int64}}(undef, 0)
+    interpolation_function_values = Vector{Vector{Float64}}(undef, 0)
+    tol = 1.0e-06
+    if haskey(bc_params, "search tolerance") == true
+        tol = bc_params["search tolerance"]
+    end
+    side_set_node_indices = unique(side_set_node_indices)
+    for node_index in side_set_node_indices
+        point = subsim.model.reference[:, node_index]
+        node_indices, ξ, found = find_point_in_mesh(point, coupled_subsim.model, coupled_block_id, tol)
+        if found == false
+            error("Could not find subdomain ", subsim.name, " point ", point, " in subdomain ", coupled_subsim.name)
+        end
+        N = interpolate(element_type, ξ)[1]
+        push!(coupled_nodes_indices, node_indices)
+        push!(interpolation_function_values, N)
+    end
+    is_dirichlet = true
+    swap_bcs = false
+
+    opinf_model_directory = bc_params["model-directory"]
+    py""" 
+    import torch
+    def get_model(model_file):
+      return torch.load(model_file)
+    """
+    ensemble_size = bc_params["ensemble-size"]
+    model = []
+    for i in 1:ensemble_size
+      tmp =  py"get_model"(opinf_model_directory * "/BC-" * side_set_name * "-" * string(i-1) * ".pt")
+      push!(model,tmp)
+    end
+
+    basis_file = bc_params["model-directory"] * "/nn-opinf-basis-" * side_set_name * ".npz" 
+    basis = NPZ.npzread(basis_file)
+    basis = basis["basis"]
+
+    if bc_type == "OpInf Schwarz overlap"
+        SMOpInfOverlapSchwarzBC(
+            side_set_name,
+            side_set_node_indices,
+            coupled_nodes_indices,
+            interpolation_function_values,
+            coupled_subsim,
+            subsim,
+            is_dirichlet,
+            swap_bcs,
+            fom_bc,
+            model,
+            basis
+        )
+    else
+        error("Unknown boundary condition type : ", bc_type)
+    end
+end
+
+
+
+function apply_bc(model::OpInfModel, bc::SMDirichletBC)
+    model.fom_model.time = model.time
+    apply_bc(model.fom_model, bc)
+    bc_vector = zeros(0)
+    for node_index in bc.node_set_node_indices
+        disp_val = model.fom_model.current[bc.offset, node_index] - model.fom_model.reference[bc.offset, node_index]
+        push!(bc_vector, disp_val)
+    end
+
+    offset = bc.offset
+    if offset == 1
+        offset_name = "x"
+    end
+    if offset == 2
+        offset_name = "y"
+    end
+    if offset == 3
+        offset_name = "z"
+    end
+
+    op_name = "B_" * bc.node_set_name * "-" * offset_name
+    bc_operator = model.opinf_rom[op_name]
+    # SM Dirichlet BC are only defined on a single x,y,z
+    return model.reduced_boundary_forcing[:] += bc_operator[1, :, :] * bc_vector
+end
+
+function apply_bc(model::NeuralNetworkOpInfRom, bc::SMOpInfDirichletBC)
+    model.fom_model.time = model.time
+    apply_bc(model.fom_model,bc.fom_bc)
+    bc_vector = zeros(0)
+    for node_index ∈ bc.fom_bc.node_set_node_indices
+        dof_index = 3 * (node_index - 1) + bc.fom_bc.offset
+        disp_val = model.fom_model.current[bc.fom_bc.offset,node_index] - model.fom_model.reference[bc.fom_bc.offset, node_index]
+        push!(bc_vector,disp_val)
+    end
+
+    py"""
+    import numpy as np
+    def setup_inputs(x):
+        xi = np.zeros((1,x.size))
+        xi[0] = x
+        inputs = torch.tensor(xi)
+        return inputs
+    """
+
+    reduced_bc_vector = bc.basis[1,:,:]' * bc_vector 
+    model_inputs = py"setup_inputs"(reduced_bc_vector)
+    ensemble_size = size(bc.nn_model)[1]
+    for i in 1:ensemble_size
+      reduced_forcing = bc.nn_model[i].forward(model_inputs)
+      reduced_forcing = reduced_forcing.detach().numpy()[1,:]
+      model.reduced_boundary_forcing[:] += reduced_forcing 
+    end
+    model.reduced_boundary_forcing[:] = model.reduced_boundary_forcing[:] ./ ensemble_size
+end
+
+function apply_bc_detail(model::OpInfModel, bc::CouplingSchwarzBoundaryCondition)
+    if (typeof(bc.coupled_subsim.model) == SolidMechanics)
+        ## Apply BC to the FOM vector
+        apply_bc_detail(model.fom_model, bc)
+
+        # populate our own BC vector
+        bc_vector = zeros(3, length(bc.side_set_node_indices))
+        for i in 1:length(bc.side_set_node_indices)
+            node_index = bc.side_set_node_indices[i]
+            bc_vector[:, i] = model.fom_model.current[:, node_index] - model.fom_model.reference[:, node_index]
+        end
+        op_name = "B_" * bc.side_set_name
+        bc_operator = model.opinf_rom[op_name]
+        for i in 1:3
+            model.reduced_boundary_forcing[:] += bc_operator[i, :, :] * bc_vector[i, :]
+        end
+    else
+        throw("ROM-ROM coupling not supported yet")
+    end
+end
+
+
+function apply_bc_detail(model::NeuralNetworkOpInfRom, bc::CouplingSchwarzBoundaryCondition)
+    if (typeof(bc.coupled_subsim.model) == SolidMechanics)
+        ## Apply BC to the FOM vector
+        apply_bc_detail(model.fom_model, bc.fom_bc)
+
+        # populate our own BC vector
+        bc_vector = zeros(3, length(bc.fom_bc.side_set_node_indices))
+        for i in 1:length(bc.fom_bc.side_set_node_indices)
+            node_index = bc.fom_bc.side_set_node_indices[i]
+            bc_vector[:, i] = model.fom_model.current[:, node_index] - model.fom_model.reference[:, node_index]
+        end  
+
+        py"""
+        import numpy as np
+        def setup_inputs(x):
+            xi = np.zeros((1,x.size))
+            xi[0] = x
+            inputs = torch.tensor(xi)
+            return inputs
+        """
+
+        reduced_bc_vector = zeros(size(bc.basis)[3])
+        for i in 1:3
+            reduced_bc_vector[:] += bc.basis[i,:,:]' * bc_vector[i,:] 
+        end  
+        model_inputs = py"setup_inputs"(reduced_bc_vector)
+        ensemble_size = size(bc.nn_model)[1]
+        local_reduced_forcing = zeros(size(model.reduced_boundary_forcing)[1])
+        for i in 1:ensemble_size
+          reduced_forcing = bc.nn_model[i].forward(model_inputs)
+          reduced_forcing = reduced_forcing.detach().numpy()[1,:]
+          local_reduced_forcing[:] += reduced_forcing 
+        end
+        local_reduced_forcing[:] = local_reduced_forcing[:] ./ ensemble_size
+        model.reduced_boundary_forcing[:] += local_reduced_forcing
+
+    else 
+        throw("ROM-ROM coupling not supported yet")
+    end  
+end
+
+function apply_bc(model::RomModel, bc::SchwarzBoundaryCondition)
+    global_sim = bc.coupled_subsim.params["global_simulation"]
+    controller = global_sim.controller
+    if bc isa SMContactSchwarzBC && controller.active_contact == false
+        return nothing
+    end
+    empty_history = length(controller.time_hist) == 0
+    same_step = controller.same_step
+    if empty_history == true
+        apply_bc_detail(model, bc)
+        return nothing
+    end
+    # Save solution of coupled simulation
+    saved_disp = bc.coupled_subsim.integrator.displacement
+    saved_velo = bc.coupled_subsim.integrator.velocity
+    saved_acce = bc.coupled_subsim.integrator.acceleration
+    saved_∂Ω_f = bc.coupled_subsim.model.internal_force
+    time = model.time
+    coupled_name = bc.coupled_subsim.name
+    coupled_index = global_sim.subsim_name_index_map[coupled_name]
+    time_hist = controller.time_hist[coupled_index]
+    disp_hist = controller.disp_hist[coupled_index]
+    velo_hist = controller.velo_hist[coupled_index]
+    acce_hist = controller.acce_hist[coupled_index]
+    ∂Ω_f_hist = controller.∂Ω_f_hist[coupled_index]
+    interp_disp = same_step == true ? disp_hist[end] : interpolate(time_hist, disp_hist, time)
+    interp_velo = same_step == true ? velo_hist[end] : interpolate(time_hist, velo_hist, time)
+    interp_acce = same_step == true ? acce_hist[end] : interpolate(time_hist, acce_hist, time)
+    interp_∂Ω_f = same_step == true ? ∂Ω_f_hist[end] : interpolate(time_hist, ∂Ω_f_hist, time)
+    if bc.coupled_subsim.model isa SolidMechanics
+        bc.coupled_subsim.model.internal_force = interp_∂Ω_f
+    elseif bc.coupled_subsim.model isa RomModel
+        bc.coupled_subsim.model.fom_model.internal_force = interp_∂Ω_f
+    end
+
+    if bc isa SMContactSchwarzBC || bc isa SMNonOverlapSchwarzBC
+        relaxation_parameter = controller.relaxation_parameter
+        schwarz_iteration = controller.iteration_number
+        if schwarz_iteration == 1
+            lambda_dispᵖʳᵉᵛ = zeros(length(interp_disp))
+            lambda_veloᵖʳᵉᵛ = zeros(length(interp_velo))
+            lambda_acceᵖʳᵉᵛ = zeros(length(interp_acce))
+        else
+            lambda_dispᵖʳᵉᵛ = controller.lambda_disp[coupled_index]
+            lambda_veloᵖʳᵉᵛ = controller.lambda_velo[coupled_index]
+            lambda_acceᵖʳᵉᵛ = controller.lambda_acce[coupled_index]
+        end
+        bc.coupled_subsim.integrator.displacement =
+            controller.lambda_disp[coupled_index] =
+                relaxation_parameter * interp_disp + (1 - relaxation_parameter) * lambda_dispᵖʳᵉᵛ
+        bc.coupled_subsim.integrator.velocity =
+            controller.lambda_velo[coupled_index] =
+                relaxation_parameter * interp_velo + (1 - relaxation_parameter) * lambda_veloᵖʳᵉᵛ
+        bc.coupled_subsim.integrator.acceleration =
+            controller.lambda_acce[coupled_index] =
+                relaxation_parameter * interp_acce + (1 - relaxation_parameter) * lambda_acceᵖʳᵉᵛ
+    else
+        bc.coupled_subsim.integrator.displacement = interp_disp
+        bc.coupled_subsim.integrator.velocity = interp_velo
+        bc.coupled_subsim.integrator.acceleration = interp_acce
+    end
+    # Copies from integrator to model
+    copy_solution_source_targets(bc.coupled_subsim.integrator, bc.coupled_subsim.solver, bc.coupled_subsim.model)
+    apply_bc_detail(model, bc)
+    bc.coupled_subsim.integrator.displacement = saved_disp
+    bc.coupled_subsim.integrator.velocity = saved_velo
+    bc.coupled_subsim.integrator.acceleration = saved_acce
+    bc.coupled_subsim.model.internal_force = saved_∂Ω_f
+    # Copy from integrator to model
+    return copy_solution_source_targets(bc.coupled_subsim.integrator, bc.coupled_subsim.solver, bc.coupled_subsim.model)
+end
+
+function apply_bcs(model::RomModel)
+    model.reduced_boundary_forcing[:] .= 0.0
+    for boundary_condition in model.boundary_conditions
+        apply_bc(model, boundary_condition)
+    end
+end
+
+
+function apply_ics(params::Parameters, model::RomModel)
+    apply_ics(params, model.fom_model)
+
+    if haskey(params, "initial conditions") == false
+        return nothing
+    end
+    n_var, n_node, n_mode = size(model.basis)
+    n_var_fom, n_node_fom = size(model.fom_model.current)
+
+    # Make sure basis is the right size
+    if n_var != n_var_fom || n_node != n_node_fom
+        throw("Basis is wrong size")
+    end
+
+    # project onto basis
+    for k in 1:n_mode
+        model.reduced_state[k] = 0.0
+        model.reduced_velocity[k] = 0.0
+        for j in 1:n_node
+            for n in 1:n_var
+                model.reduced_state[k] +=
+                    model.basis[n, j, k] *
+                    (model.fom_model.current[n, j] - model.fom_model.reference[n, j])
+                model.reduced_velocity[k] +=
+                    model.basis[n, j, k] *
+                    (model.fom_model.velocity[n, j])
+            end
+        end
+    end
+end
+
+
+function find_point_in_mesh(point::Vector{Float64}, model::RomModel, blk_id::Int, tol::Float64)
+    node_indices, ξ, found = find_point_in_mesh(point, model.fom_model, blk_id, tol)
+    return node_indices, ξ, found
+end
+
diff --git a/src/opinf/opinf_ics_bcs_types.jl b/src/opinf/opinf_ics_bcs_types.jl
new file mode 100644
index 00000000..f9deea1e
--- /dev/null
+++ b/src/opinf/opinf_ics_bcs_types.jl
@@ -0,0 +1,47 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+
+abstract type BoundaryCondition end
+abstract type SchwarzBoundaryCondition <: BoundaryCondition end
+abstract type RegularBoundaryCondition <: BoundaryCondition end
+abstract type ContactSchwarzBoundaryCondition <: SchwarzBoundaryCondition end
+abstract type CouplingSchwarzBoundaryCondition <: SchwarzBoundaryCondition end
+abstract type OverlapSchwarzBoundaryCondition <: CouplingSchwarzBoundaryCondition end
+abstract type NonOverlapSchwarzBoundaryCondition <: CouplingSchwarzBoundaryCondition end
+abstract type InitialCondition end
+
+using Exodus
+using Symbolics
+
+mutable struct SMOpInfDirichletBC <: RegularBoundaryCondition 
+    node_set_name::String
+    offset::Int64
+    node_set_id::Int64
+    node_set_node_indices::Vector{Int64}
+    disp_num::Num
+    velo_num::Num
+    acce_num::Num
+    fom_bc::SMDirichletBC
+    nn_model::Any
+    basis::Array{Float64}
+end
+
+
+mutable struct SMOpInfOverlapSchwarzBC <: OverlapSchwarzBoundaryCondition
+    side_set_name::String
+    side_set_node_indices::Vector{Int64}
+    coupled_nodes_indices::Vector{Vector{Int64}}
+    interpolation_function_values::Vector{Vector{Float64}}
+    coupled_subsim::Simulation
+    subsim::Simulation
+    is_dirichlet::Bool
+    swap_bcs::Bool
+    fom_bc::SMOverlapSchwarzBC
+    nn_model::Any
+    basis::Array{Float64}
+end
+
+
diff --git a/src/opinf/opinf_model.jl b/src/opinf/opinf_model.jl
new file mode 100644
index 00000000..8f17d18b
--- /dev/null
+++ b/src/opinf/opinf_model.jl
@@ -0,0 +1,129 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+
+using Base.Threads: @threads, threadid, nthreads
+using NPZ
+using PyCall
+
+function LinearOpInfRom(params::Parameters)
+    params["mesh smoothing"] = false
+    fom_model = SolidMechanics(params)
+    reference = fom_model.reference
+    opinf_model_file = params["model"]["model-file"]
+    opinf_model = NPZ.npzread(opinf_model_file)
+    basis = opinf_model["basis"]
+    _, _, reduced_dim = size(basis)
+    num_dofs = reduced_dim
+    time = 0.0
+    failed = false
+    null_vec = zeros(num_dofs)
+
+    reduced_state = zeros(num_dofs)
+    reduced_velocity = zeros(num_dofs)
+    reduced_boundary_forcing = zeros(num_dofs)
+    free_dofs = trues(num_dofs)
+    boundary_conditions = Vector{BoundaryCondition}()
+    return LinearOpInfRom(
+        opinf_model,
+        basis,
+        reduced_state,
+        reduced_velocity,
+        reduced_boundary_forcing,
+        null_vec,
+        free_dofs,
+        boundary_conditions,
+        time,
+        failed,
+        fom_model,
+        reference,
+        false,
+    )
+end
+
+function QuadraticOpInfRom(params::Parameters)
+    params["mesh smoothing"] = false
+    fom_model = SolidMechanics(params)
+    reference = fom_model.reference
+    opinf_model_file = params["model"]["model-file"]
+    opinf_model = NPZ.npzread(opinf_model_file)
+    basis = opinf_model["basis"]
+    _, _, reduced_dim = size(basis)
+    num_dofs = reduced_dim
+    time = 0.0
+    failed = false
+    null_vec = zeros(num_dofs)
+
+    reduced_state = zeros(num_dofs)
+    reduced_velocity = zeros(num_dofs)
+    reduced_boundary_forcing = zeros(num_dofs)
+    free_dofs = trues(num_dofs)
+    boundary_conditions = Vector{BoundaryCondition}()
+    return QuadraticOpInfRom(
+        opinf_model,
+        basis,
+        reduced_state,
+        reduced_velocity,
+        reduced_boundary_forcing,
+        null_vec,
+        free_dofs,
+        boundary_conditions,
+        time,
+        failed,
+        fom_model,
+        reference,
+        false,
+    )
+end
+
+function NeuralNetworkOpInfRom(params::Dict{String,Any})
+    params["mesh smoothing"] = false
+    fom_model = SolidMechanics(params)
+    reference = fom_model.reference
+    opinf_model_directory = params["model"]["model-directory"]
+    basis_file = opinf_model_directory * "/nn-opinf-basis.npz"
+    basis = NPZ.npzread(basis_file)
+    basis = basis["basis"]
+    py""" 
+    import torch
+    def get_model(model_file):
+      return torch.load(model_file)
+    """
+    ensemble_size = params["model"]["ensemble-size"]
+    model = []
+    for i in 1:ensemble_size
+      tmp =  py"get_model"(opinf_model_directory * "/stiffness-" * string(i-1) * ".pt")
+      push!(model,tmp)
+    end
+    num_dofs_per_node,num_nodes_basis,reduced_dim = size(basis)
+    num_dofs = reduced_dim
+
+    time = 0.0
+    failed = false
+    null_vec = zeros(num_dofs)
+
+    reduced_state = zeros(num_dofs)
+    reduced_velocity = zeros(num_dofs)
+    reduced_boundary_forcing = zeros(num_dofs)
+    free_dofs = trues(num_dofs)
+    boundary_conditions = Vector{BoundaryCondition}()
+    NeuralNetworkOpInfRom(
+        model,
+        basis,
+        reduced_state,
+        reduced_velocity,
+        reduced_boundary_forcing,
+        null_vec,
+        free_dofs,
+        boundary_conditions,
+        time,
+        failed,
+        fom_model,
+        reference,
+        false 
+    )
+end
+
+
diff --git a/src/opinf/opinf_model_types.jl b/src/opinf/opinf_model_types.jl
new file mode 100644
index 00000000..29963da8
--- /dev/null
+++ b/src/opinf/opinf_model_types.jl
@@ -0,0 +1,60 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+
+abstract type RomModel <: Model end
+abstract type OpInfModel <: RomModel end
+
+mutable struct QuadraticOpInfRom <: OpInfModel
+    opinf_rom::Dict{Any,Any}
+    basis::Array{Float64}
+    reduced_state::Vector{Float64}
+    reduced_velocity::Vector{Float64}
+    reduced_boundary_forcing::Vector{Float64}
+    #internal_force not used, but include to ease interfacing in Schwarz
+    internal_force::Vector{Float64}
+    free_dofs::BitVector
+    boundary_conditions::Vector{BoundaryCondition}
+    time::Float64
+    failed::Bool
+    fom_model::SolidMechanics
+    reference::Matrix{Float64}
+    inclined_support::Bool
+end
+
+
+mutable struct LinearOpInfRom <: OpInfModel
+    opinf_rom::Dict{Any,Any}
+    basis::Array{Float64}
+    reduced_state::Vector{Float64}
+    reduced_velocity::Vector{Float64}
+    reduced_boundary_forcing::Vector{Float64}
+    #internal_force not used, but include to ease interfacing in Schwarz
+    internal_force::Vector{Float64}
+    free_dofs::BitVector
+    boundary_conditions::Vector{BoundaryCondition}
+    time::Float64
+    failed::Bool
+    fom_model::SolidMechanics
+    reference::Matrix{Float64}
+    inclined_support::Bool
+end
+
+mutable struct NeuralNetworkOpInfRom <: OpInfModel
+    nn_model::Any
+    basis::Array{Float64}
+    reduced_state::Vector{Float64}
+    reduced_velocity::Vector{Float64}
+    reduced_boundary_forcing::Vector{Float64}
+    #internal_force not used, but include to ease interfacing in Schwarz
+    internal_force::Vector{Float64}
+    free_dofs::BitVector
+    boundary_conditions::Vector{BoundaryCondition}
+    time::Float64
+    failed::Bool
+    fom_model::SolidMechanics
+    reference::Matrix{Float64}
+    inclined_support::Bool
+end
diff --git a/src/simulation.jl b/src/simulation.jl
index ae8d91c6..df35bd05 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -14,6 +14,9 @@ include("model.jl")
 include("time_integrator.jl")
 include("solver.jl")
 include("schwarz.jl")
+# For OpInf models
+include("opinf/opinf_model.jl")
+
 
 function create_simulation(input_file::String)
     norma_log(0, :setup, "Reading from " * input_file)
diff --git a/src/simulation_types.jl b/src/simulation_types.jl
index 9669db85..bdcf8f04 100644
--- a/src/simulation_types.jl
+++ b/src/simulation_types.jl
@@ -12,6 +12,7 @@ include("ics_bcs_types.jl")
 include("model_types.jl")
 include("time_integrator_types.jl")
 include("solver_types.jl")
+include("opinf/opinf_model_types.jl")
 
 abstract type SingleController end
 abstract type MultiDomainController end
@@ -89,3 +90,4 @@ mutable struct MultiDomainSimulation <: Simulation
     subsim_name_index_map::Dict{String,Int64}
     failed::Bool
 end
+
diff --git a/src/solver.jl b/src/solver.jl
index bde7568b..fc691106 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -408,15 +408,6 @@ function evaluate(integrator::Newmark, solver::HessianMinimizer, model::LinearOp
     beta = integrator.β
     gamma = integrator.γ
     dt = integrator.time_step
-
-    #Ax* = b
-    # Put in residual format
-    #e = [x* - x] -> x* = x + e
-    #Ax + Ae = b
-    #Ax - b = -Ae
-    #Ae = r, r = b - Ax 
-    ##M uddot + Ku = f
-
     num_dof = length(model.free_dofs)
     I = Matrix{Float64}(LinearAlgebra.I, num_dof, num_dof)
     LHS = I / (dt * dt * beta) + Matrix{Float64}(model.opinf_rom["K"])
@@ -428,6 +419,44 @@ function evaluate(integrator::Newmark, solver::HessianMinimizer, model::LinearOp
     return nothing
 end
 
+### Move to model?  
+function evaluate(integrator::Newmark, solver::HessianMinimizer, model::NeuralNetworkOpInfRom)
+    beta  = integrator.β
+    gamma = integrator.γ
+    dt = integrator.time_step
+    num_dof, = size(model.free_dofs)
+    I = Matrix{Float64}(LinearAlgebra.I, num_dof,num_dof)
+    py"""
+    import numpy as np
+    def setup_inputs(x):
+        xi = np.zeros((1,x.size))
+        xi[0] = x
+        inputs = torch.tensor(xi)
+        return inputs
+    """
+    model_inputs = py"setup_inputs"(solver.solution)
+    ensemble_size = size(model.nn_model)[1] 
+    stiffness = zeros( num_dof,num_dof )
+    #Kx,K = model.nn_model[1].forward(model_inputs,return_stiffness=true)
+    #K = K.detach().numpy()[1,:,:]
+    for i in 1:ensemble_size
+      Kxt,Kt = model.nn_model[i].forward(model_inputs,return_stiffness=true)
+      #Kx += Kxt
+      Kt = Kt.detach().numpy()[1,:,:]
+      stiffness += Kt
+
+    end
+    stiffness = stiffness./ensemble_size
+    LHS = I / (dt*dt*beta) - stiffness 
+    RHS = model.reduced_boundary_forcing + 1.0/(dt*dt*beta).*integrator.disp_pre
+
+    residual = RHS - LHS * solver.solution 
+    solver.hessian[:,:] = LHS
+    solver.gradient[:] = -residual 
+end
+
+
+
 function evaluate(integrator::QuasiStatic, solver::HessianMinimizer, model::SolidMechanics)
     evaluate(model, integrator, solver)
     if model.failed == true
diff --git a/test/nnopinf-schwarz-overlap-cuboid-hex8.jl b/test/nnopinf-schwarz-overlap-cuboid-hex8.jl
new file mode 100644
index 00000000..bd167006
--- /dev/null
+++ b/test/nnopinf-schwarz-overlap-cuboid-hex8.jl
@@ -0,0 +1,29 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+@testset "Neural network Opinf Schwarz Overlap Static Cuboid Hex8 Same Step" begin
+    cp("../examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-1-for-test.yaml", "cuboid-1-for-test.yaml"; force=true)
+    cp("../examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboid-2-for-test.yaml", "cuboid-2-for-test.yaml"; force=true)
+    cp("../examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/cuboids-for-test.yaml", "cuboids-for-test.yaml"; force=true)
+    cp("../examples/ahead/overlap/cuboid/dynamic-opinf-nn-rom/neural-network-model", "neural-network-model"; force=true)
+    cp("../examples/ahead/overlap/cuboid/dynamic-opinf-rom/cuboid-1.g", "cuboid-1.g"; force=true)
+    cp("../examples/ahead/overlap/cuboid/dynamic-opinf-rom/cuboid-2.g", "cuboid-2.g"; force=true)
+    sim = Norma.run("cuboids-for-test.yaml")
+    subsims = sim.subsims
+    model_fine = subsims[1].model
+    model_coarse = subsims[2].model
+    rm("cuboids-for-test.yaml")
+    rm("cuboid-1-for-test.yaml")
+    rm("cuboid-2-for-test.yaml")
+    rm("cuboid-1.g")
+    rm("cuboid-2.g")
+    rm("cuboid-1.e")
+    rm("cuboid-2.e")
+    for file in readdir("neural-network-model")
+      rm(joinpath("neural-network-model", file); force=true)  # Remove each file/directory
+    end
+
+    rm("neural-network-model"; force=true)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7daccf63..103de30a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -12,7 +12,7 @@ using Test
 include("../src/Norma.jl")
 include("helpers.jl")
 
-const all_test_files = [
+all_test_files = [
     "minitensor.jl",
     "interpolation.jl",
     "constitutive.jl",
@@ -35,10 +35,20 @@ const all_test_files = [
     "opinf-schwarz-overlap-cuboid-hex8.jl",
     "quadratic-opinf-schwarz-overlap-cuboid-hex8.jl",
     "adaptive-time-stepping.jl",
-    # Must go last for now due to FPE trapping
-    "utils.jl",
 ]
 
+nn_opinf_test_files = [
+    "nnopinf-schwarz-overlap-cuboid-hex8.jl",
+]
+
+if "run-nn-opinf-tests" in ARGS
+   append!(all_test_files,nn_opinf_test_files)
+end
+
+# Must go last for now due to FPE trapping
+"utils.jl",
+push!(all_test_files,"utils.jl")
+
 const indexed_test_files = collect(enumerate(all_test_files))
 
 # Display test options
@@ -48,8 +58,9 @@ for (i, file) in indexed_test_files
 end
 
 # Parse command-line arguments as indices (if any)
-selected_test_indices = parse.(Int, ARGS)
 
+selected_test_indices = parse.(Int, filter(arg-> arg!= "run-nn-opinf-tests",ARGS))
+       
 # Determine which tests to run
 test_files_to_run = isempty(selected_test_indices) ? indexed_test_files :
     filter(t -> t[1] in selected_test_indices, indexed_test_files)