diff --git a/Source/C/DensityMatrixSolvers_c.h b/Source/C/DensityMatrixSolvers_c.h
index 0df5ebc1..2ab58d43 100644
--- a/Source/C/DensityMatrixSolvers_c.h
+++ b/Source/C/DensityMatrixSolvers_c.h
@@ -33,4 +33,7 @@ void DenseDensity_wrp(const int *ih_Hamiltonian,
                       const int *ih_solver_parameters);
 void EnergyDensityMatrix_wrp(const int *ih_Hamiltonian, const int *ih_Density,
                              int *ih_EnergyDensity, const double *threshold);
+void McWeenyStep_wrp(const int *ih_D, int *ih_DOut, const double *threshold);
+void McWeenyStepS_wrp(const int *ih_D, int *ih_DOut, const int *ih_S,
+                      const double *threshold);
 #endif
diff --git a/Source/CPlusPlus/DensityMatrixSolvers.cc b/Source/CPlusPlus/DensityMatrixSolvers.cc
index 29222494..9aff0c1f 100644
--- a/Source/CPlusPlus/DensityMatrixSolvers.cc
+++ b/Source/CPlusPlus/DensityMatrixSolvers.cc
@@ -80,4 +80,16 @@ void DensityMatrixSolvers::EnergyDensityMatrix(const Matrix_ps &Hamiltonian,
                           GetIH(EnergyDensity), &threshold);
 }
 
+////////////////////////////////////////////////////////////////////////////////
+void DensityMatrixSolvers::McWeenyStep(const Matrix_ps &D, Matrix_ps &DOut,
+                                       double threshold) {
+  McWeenyStep_wrp(GetIH(D), GetIH(DOut), &threshold);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void DensityMatrixSolvers::McWeenyStep(const Matrix_ps &D, const Matrix_ps &S,
+                                       Matrix_ps &DOut, double threshold) {
+  McWeenyStepS_wrp(GetIH(D), GetIH(DOut), GetIH(S), &threshold);
+}
+
 } // namespace NTPoly
diff --git a/Source/CPlusPlus/DensityMatrixSolvers.h b/Source/CPlusPlus/DensityMatrixSolvers.h
index 2d7009fa..fc84500c 100644
--- a/Source/CPlusPlus/DensityMatrixSolvers.h
+++ b/Source/CPlusPlus/DensityMatrixSolvers.h
@@ -106,6 +106,19 @@ class DensityMatrixSolvers : public SolverBase {
                                   const Matrix_ps &Density,
                                   Matrix_ps &EnergyDensity,
                                   double threshold = 0.0);
+  //! Take one McWeeny Step DOut = 3*DD - 2*DDD
+  //!\param D the density matrix.
+  //!\param DOut a more purified matrix.
+  //!\param threshold for flushing small values to zero.
+  static void McWeenyStep(const Matrix_ps &D, Matrix_ps &DOut,
+                          double threshold = 0.0);
+  //! Take one McWeeny Step DOut = 3*DSD - 2*DSDSD
+  //!\param D the density matrix.
+  //!\param S the overlap matrix.
+  //!\param DOut a more purified matrix.
+  //!\param threshold for flushing small values to zero.
+  static void McWeenyStep(const Matrix_ps &D, const Matrix_ps &S,
+                          Matrix_ps &DOut, double threshold = 0.0);
 };
 } // namespace NTPoly
 #endif
diff --git a/Source/Fortran/DensityMatrixSolversModule.F90 b/Source/Fortran/DensityMatrixSolversModule.F90
index c1be31f9..a118527f 100644
--- a/Source/Fortran/DensityMatrixSolversModule.F90
+++ b/Source/Fortran/DensityMatrixSolversModule.F90
@@ -28,6 +28,7 @@ MODULE DensityMatrixSolversModule
   PUBLIC :: DenseDensity
   PUBLIC :: ScaleAndFold
   PUBLIC :: EnergyDensityMatrix
+  PUBLIC :: McWeenyStep
 CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Compute the density matrix from a Hamiltonian using the PM method.
   !> Based on the PM algorithm presented in \cite palser1998canonical
@@ -1207,5 +1208,49 @@ SUBROUTINE EnergyDensityMatrix(H, D, ED, threshold_in)
     CALL SimilarityTransform(H, D, D, ED, threshold_in=threshold)
 
   END SUBROUTINE EnergyDensityMatrix
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Take one McWeeny Step DOut = 3*DSD - 2*DSDSD
+  SUBROUTINE McWeenyStep(D, DOut, S_in, threshold_in)
+    !> The density matrix.
+    TYPE(Matrix_ps), INTENT(IN) :: D
+    !> The resulting purified matrix.
+    TYPE(Matrix_ps), INTENT(INOUT) :: DOut
+    !> The overlap matrix (optional)
+    TYPE(Matrix_ps), INTENT(IN), OPTIONAL :: S_in
+    !> Threshold for flushing small values (default = 0).
+    REAL(NTREAL), INTENT(IN), OPTIONAL :: threshold_in
+    !! Handling Optional Parameters
+    REAL(NTREAL) :: threshold
+    !! Local Variables
+    TYPE(Matrix_ps) :: DS, DSD
+    TYPE(MatrixMemoryPool_p) :: pool
+
+    !! Optional Parameters
+    IF (PRESENT(threshold_in)) THEN
+       threshold = threshold_in
+    ELSE
+       threshold = 0.0_NTREAL
+    END IF
+
+    !! Form the matrix DS
+    IF (PRESENT(S_in)) THEN
+       CALL MatrixMultiply(D, S_in, DS, &
+            & threshold_in=threshold, memory_pool_in=pool)
+    ELSE
+       CALL CopyMatrix(D, DS)
+    END IF
+
+    !! Compute
+    CALL MatrixMultiply(DS, D, DSD, & 
+         & threshold_in=threshold, memory_pool_in=pool)
+    CALL MatrixMultiply(DS, DSD, DOut, alpha_in=-2.0_NTREAL, &
+         & threshold_in=threshold, memory_pool_in=pool)
+    CALL IncrementMatrix(DSD, DOut, alpha_in=3.0_NTREAL)
+
+    !! Cleanup
+    CALL DestructMatrix(DS)
+    CALL DestructMatrix(DSD)
+    CALL DestructMatrixMemoryPool(pool)
+  END SUBROUTINE McWeenyStep
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 END MODULE DensityMatrixSolversModule
diff --git a/Source/Wrapper/DensityMatrixSolversModule_wrp.F90 b/Source/Wrapper/DensityMatrixSolversModule_wrp.F90
index 54e4cdd5..ed9b478f 100644
--- a/Source/Wrapper/DensityMatrixSolversModule_wrp.F90
+++ b/Source/Wrapper/DensityMatrixSolversModule_wrp.F90
@@ -17,6 +17,7 @@ MODULE DensityMatrixSolversModule_wrp
   PUBLIC :: DenseDensity_wrp
   PUBLIC :: ScaleAndFold_wrp
   PUBLIC :: EnergyDensityMatrix_wrp
+  PUBLIC :: McWeenyStep_wrp
   ! PUBLIC :: HPCPPlus_wrp
 CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Compute the density matrix from a Hamiltonian using the PM method.
@@ -196,5 +197,38 @@ SUBROUTINE EnergyDensityMatrix_wrp(ih_Hamiltonian, ih_Density, &
     CALL EnergyDensityMatrix(h_Hamiltonian%DATA, h_Density%DATA, &
          & h_EnergyDensity%DATA, threshold)
   END SUBROUTINE EnergyDensityMatrix_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Take one McWeeny Step DOut = 3*DD - 2*DDD
+  SUBROUTINE McWeenyStep_wrp(ih_D, ih_DOut, threshold) &
+       & BIND(c,name="McWeenyStep_wrp")
+    INTEGER(kind=c_int), INTENT(IN) :: ih_D(SIZE_wrp)
+    INTEGER(kind=c_int), INTENT(INOUT) :: ih_DOut(SIZE_wrp)
+    REAL(NTREAL), INTENT(IN) :: threshold
+    TYPE(Matrix_ps_wrp) :: h_D
+    TYPE(Matrix_ps_wrp) :: h_DOut
+
+    h_D = TRANSFER(ih_D,h_D)
+    h_DOut = TRANSFER(ih_DOut,h_DOut)
+
+    CALL McWeenyStep(h_D%DATA, h_Dout%DATA, threshold_in=threshold)
+  END SUBROUTINE McWeenyStep_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Take one McWeeny Step DOut = 3*DSD - 2*DSDSD
+  SUBROUTINE McWeenyStepS_wrp(ih_D, ih_DOut, ih_S, threshold) &
+       & BIND(c,name="McWeenyStepS_wrp")
+    INTEGER(kind=c_int), INTENT(IN) :: ih_D(SIZE_wrp)
+    INTEGER(kind=c_int), INTENT(INOUT) :: ih_DOut(SIZE_wrp)
+    INTEGER(kind=c_int), INTENT(IN) :: ih_S(SIZE_wrp)
+    REAL(NTREAL), INTENT(IN) :: threshold
+    TYPE(Matrix_ps_wrp) :: h_D, h_S
+    TYPE(Matrix_ps_wrp) :: h_DOut
+
+    h_D = TRANSFER(ih_D,h_D)
+    h_S = TRANSFER(ih_S,h_S)
+    h_DOut = TRANSFER(ih_DOut,h_DOut)
+
+    CALL McWeenyStep(h_D%DATA, h_Dout%DATA, S_in=h_S%DATA, &
+         & threshold_in=threshold)
+  END SUBROUTINE McWeenyStepS_wrp
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 END MODULE DensityMatrixSolversModule_wrp
diff --git a/UnitTests/test_chemistry.py b/UnitTests/test_chemistry.py
index 3927617f..b8408eb1 100644
--- a/UnitTests/test_chemistry.py
+++ b/UnitTests/test_chemistry.py
@@ -287,6 +287,49 @@ def test_foe_low(self):
         '''Test the fermi operator expansion at a low temperature.'''
         self.basic_solver(nt.FermiOperator.ComputeDenseFOE, temp=50)
 
+    def test_mcweeny_step(self):
+        from scipy.io import mmread
+        from helpers import result_file, THRESHOLD
+        from scipy.sparse.linalg import norm
+
+        # Reference Solution
+        dmat = mmread(self.density)
+        smat = mmread(self.overlap)
+        result = 3 * dmat.dot(dmat) - 2 * dmat.dot(dmat).dot(dmat)
+        result_s = 3 * dmat.dot(smat).dot(dmat) - \
+            2 * dmat.dot(smat).dot(dmat).dot(smat).dot(dmat)
+
+        # NTPoly
+        d_matrix = nt.Matrix_ps(self.density)
+        dout_matrix = nt.Matrix_ps(d_matrix.GetActualDimension())
+        nt.DensityMatrixSolvers.McWeenyStep(d_matrix, dout_matrix)
+        dout_matrix.WriteToMatrixMarket(result_file)
+        comm.barrier()
+
+        # Compare
+        normval = 0
+        if (self.my_rank == 0):
+            ResultMat = mmread(result_file)
+            normval = abs(norm(result - ResultMat))
+        global_norm = comm.bcast(normval, root=0)
+        self.assertLessEqual(global_norm, THRESHOLD)
+        comm.barrier()
+
+        # Overlap Version
+        s_matrix = nt.Matrix_ps(self.overlap)
+        nt.DensityMatrixSolvers.McWeenyStep(d_matrix, s_matrix, dout_matrix)
+        dout_matrix.WriteToMatrixMarket(result_file)
+        comm.barrier()
+
+        # Compare
+        normval = 0
+        if (self.my_rank == 0):
+            ResultMat = mmread(result_file)
+            normval = abs(norm(result_s - ResultMat))
+        global_norm = comm.bcast(normval, root=0)
+        self.assertLessEqual(global_norm, THRESHOLD)
+        comm.barrier()
+
     def test_energy_density(self):
         '''Test the routines to compute the weighted-energy density matrix.'''
         from helpers import THRESHOLD, result_file