diff --git a/src/qs_scf.F b/src/qs_scf.F
index edcf070e07..a8141370ce 100644
--- a/src/qs_scf.F
+++ b/src/qs_scf.F
@@ -480,8 +480,14 @@ CONTAINS
                                                diis_step,use_jacobi,error)
             END IF
           CASE(special_diag_method_nr)
-             CALL qs_scf_loop_do_special_diag(scf_env,mos,rho,matrix_ks,scf_control,scf_section, &
-                                          diis_step,error)
+            IF (dft_control%roks) THEN
+              CALL qs_scf_loop_do_roks_diag(scf_env,mos,matrix_ks,matrix_s,&
+                                            scf_control,scf_section,diis_step,&
+                                            has_unit_metric,error)
+            ELSE
+              CALL qs_scf_loop_do_special_diag(scf_env,mos,rho,matrix_ks,scf_control,scf_section, &
+                                               diis_step,error)
+            END IF
           CASE(ot_method_nr) 
              CALL qs_scf_loop_do_ot(scf_env,mos,rho,qs_env%mo_derivs,energy%total, & 
                                     matrix_s, energy_only, has_unit_metric,error)
@@ -2223,17 +2229,21 @@ END SUBROUTINE qs_scf_print_summary
 
     CALL cp_fm_syevd(ksa,work,eiga)
 
-    ! Use old MO set for back-transformation if level-shifting was applied
-
-    IF (level_shift_loc /= 0.0_dp) CALL cp_fm_to_fm(moa,ortho)
-
     ! Back-transformation of the orthonormal eigenvectors if needed
 
-    IF ((level_shift_loc /= 0.0_dp).OR.(.NOT.orthogonal_basis)) THEN
-      IF (scf_control%use_cholesky) THEN
-        CALL cp_fm_cholesky_restore(work,nmo,ortho,moa,"SOLVE")
+    IF (level_shift_loc /= 0.0_dp) THEN
+      ! Use old MO set for back-transformation if level-shifting was applied
+      CALL cp_fm_to_fm(moa,ortho)
+      CALL cp_fm_gemm("N","N",nao,nmo,nao,1.0_dp,ortho,work,0.0_dp,moa)
+    ELSE
+      IF (orthogonal_basis) THEN
+        CALL cp_fm_to_fm(work,moa)
       ELSE
-        CALL cp_fm_gemm("N","N",nao,nmo,nao,1.0_dp,ortho,work,0.0_dp,moa)
+        IF (scf_control%use_cholesky) THEN
+          CALL cp_fm_cholesky_restore(work,nmo,ortho,moa,"SOLVE")
+        ELSE
+          CALL cp_fm_gemm("N","N",nao,nmo,nao,1.0_dp,ortho,work,0.0_dp,moa)
+        END IF
       END IF
     END IF
 
diff --git a/tests/SE/regtest/NO2-ROKS.inp b/tests/SE/regtest/NO2-ROKS.inp
new file mode 100644
index 0000000000..75d5d8a933
--- /dev/null
+++ b/tests/SE/regtest/NO2-ROKS.inp
@@ -0,0 +1,55 @@
+&GLOBAL
+  PROGRAM Quickstep
+  PROJECT NO2
+  RUN_TYPE GEO_OPT
+  PRINT_LEVEL LOW
+&END GLOBAL
+
+&MOTION
+ &GEOOPT
+  OPTIMIZER BFGS
+  MAX_ITER  20
+  MAX_DR    0.001
+  RMS_DR    0.0005
+  MAX_FORCE 0.00015
+  RMS_FORCE 0.0001
+ &END GEOOPT
+&END MOTION
+
+&FORCE_EVAL
+  &DFT
+    ROKS
+    MULTIP 2
+    &QS
+      EXTRAPOLATION use_prev_p
+      METHOD PM3
+    &END QS
+    &SCF
+      EPS_SCF 1.0E-5
+      SCF_GUESS atomic
+    &END SCF
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 10.0 10.0 10.0
+      UNIT ANGSTROM
+      PERIODIC NONE
+    &END CELL
+    &COORD
+      N     0.000000     0.000000     0.000000
+      O     0.000000     0.000000     1.194550
+      O     0.861442     0.000000    -0.827567
+    &END COORD
+    &PRINT
+      &STRUCTURE_DATA
+        DISTANCE 1 2
+        DISTANCE 1 3
+        ANGLE 2 1 3
+      &END STRUCTURE_DATA
+      &KINDS
+        BASIS_SET
+        POTENTIAL
+      &END KINDS
+    &END PRINT
+  &END SUBSYS
+&END FORCE_EVAL
diff --git a/tests/SE/regtest/O-ROKS.inp b/tests/SE/regtest/O-ROKS.inp
new file mode 100644
index 0000000000..10e17a0a86
--- /dev/null
+++ b/tests/SE/regtest/O-ROKS.inp
@@ -0,0 +1,43 @@
+&GLOBAL
+  PROGRAM Quickstep
+  PROJECT O
+  RUN_TYPE ENERGY
+  PRINT_LEVEL LOW
+&END GLOBAL
+
+&FORCE_EVAL
+  &DFT
+    ROKS
+    MULTIP 3
+    &PRINT
+      &MO
+        EIGENVALUES
+        EIGENVECTORS
+        OCCUPATION_NUMBERS
+      &END MO
+    &END PRINT
+    &QS
+      METHOD AM1
+    &END QS
+    &SCF
+      EPS_SCF 1.0E-6
+      SCF_GUESS CORE
+    &END SCF
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 10.0 10.0 10.0
+      UNIT ANGSTROM
+      PERIODIC NONE
+    &END CELL
+    &COORD
+      O     0.000000     0.000000     0.000000
+    &END COORD
+    &PRINT
+      &KINDS
+        BASIS_SET
+        POTENTIAL
+      &END KINDS
+    &END PRINT
+  &END SUBSYS
+&END FORCE_EVAL
diff --git a/tests/SE/regtest/O2-ROKS.inp b/tests/SE/regtest/O2-ROKS.inp
new file mode 100644
index 0000000000..85cbd7dea4
--- /dev/null
+++ b/tests/SE/regtest/O2-ROKS.inp
@@ -0,0 +1,51 @@
+&GLOBAL
+  PROGRAM Quickstep
+  PROJECT O2
+  RUN_TYPE GEO_OPT
+  PRINT_LEVEL LOW
+&END GLOBAL
+
+&MOTION
+ &GEOOPT
+  OPTIMIZER BFGS
+  MAX_ITER  30
+  MAX_DR    0.001
+  RMS_DR    0.0005
+  MAX_FORCE 0.00015
+  RMS_FORCE 0.0001
+ &END GEOOPT
+&END MOTION
+
+&FORCE_EVAL
+  &DFT
+    ROKS
+    MULTIP 3
+    &QS
+      METHOD MNDO
+    &END QS
+    &SCF
+      EPS_SCF 1.0E-6
+      SCF_GUESS atomic
+    &END SCF
+  &END DFT
+  &SUBSYS
+    &CELL
+      ABC 10.0 10.0 10.0
+      UNIT ANGSTROM
+      PERIODIC NONE
+    &END CELL
+    &COORD
+      O     0.000000     0.000000    -0.600000
+      O     0.000000     0.000000     0.600000
+    &END COORD
+    &PRINT
+      &STRUCTURE_DATA
+        DISTANCE 1 2
+      &END STRUCTURE_DATA
+      &KINDS
+        BASIS_SET
+        POTENTIAL
+      &END KINDS
+    &END PRINT
+  &END SUBSYS
+&END FORCE_EVAL
diff --git a/tests/SE/regtest/TEST_FILES b/tests/SE/regtest/TEST_FILES
index 67da75f6f2..efc59f9113 100644
--- a/tests/SE/regtest/TEST_FILES
+++ b/tests/SE/regtest/TEST_FILES
@@ -18,3 +18,7 @@ st.inp 3
 # test the initial adaptation to OT
 ch4-ot.inp 3
 h2o_lsd-ot.inp 3
+# tests for high-spin ROKS
+O-ROKS.inp 3
+O2-ROKS.inp 3
+NO2-ROKS.inp 3