diff --git a/CMakeLists.txt b/CMakeLists.txt
index 753a455d..4d16166a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -4,7 +4,7 @@ cmake_policy(SET CMP0074 NEW)
 
 project (MGmol C CXX Fortran)
 
-set (CMAKE_CXX_STANDARD 11)
+set (CMAKE_CXX_STANDARD 14)
 
 # Specify the location of additional CMAKE modules
 SET(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake_modules)
diff --git a/drivers/CMakeLists.txt b/drivers/CMakeLists.txt
index 4eb1c741..36182df2 100644
--- a/drivers/CMakeLists.txt
+++ b/drivers/CMakeLists.txt
@@ -4,8 +4,14 @@ add_executable(check_input check_input.cc)
 
 add_executable(example1 example1.cc)
 
+add_executable(lbfgspp lbfgspp.cc)
+
 target_include_directories(check_input PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(example1 PRIVATE ${Boost_INCLUDE_DIRS})
+target_include_directories(lbfgspp PRIVATE ${Boost_INCLUDE_DIRS})
+target_include_directories(lbfgspp PRIVATE /home/q8j/GIT/eigen)
+target_include_directories(lbfgspp PRIVATE /home/q8j/GIT/LBFGSpp/include)
 
 target_link_libraries(check_input mgmol_src)
 target_link_libraries(example1 mgmol_src)
+target_link_libraries(lbfgspp mgmol_src)
diff --git a/drivers/lbfgspp.cc b/drivers/lbfgspp.cc
new file mode 100644
index 00000000..6ee8be54
--- /dev/null
+++ b/drivers/lbfgspp.cc
@@ -0,0 +1,249 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+
+// This driver uses the header-only LBFGS++ library
+// https://github.com/yixuan/LBFGSpp
+
+#include <Eigen/Core>
+#include <LBFGS.h>
+#include <iostream>
+
+#include "Control.h"
+#include "ExtendedGridOrbitals.h"
+#include "LocGridOrbitals.h"
+#include "MGmol.h"
+#include "MGmol_MPI.h"
+#include "MPIdata.h"
+#include "mgmol_run.h"
+
+#include <cassert>
+#include <iostream>
+#include <time.h>
+#include <vector>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+using Eigen::VectorXd;
+
+class MGmolEnergyAndForces
+{
+private:
+    MGmolInterface* mgmol_;
+    int n_;
+    std::vector<short>& anumbers_;
+
+public:
+    MGmolEnergyAndForces(
+        MGmolInterface* mgmol, const int n, std::vector<short>& anumbers)
+        : mgmol_(mgmol), n_(n), anumbers_(anumbers)
+    {
+    }
+
+    double operator()(const VectorXd& x, VectorXd& grad)
+    {
+        std::vector<double> positions(n_);
+        for (int i = 0; i < n_; i++)
+        {
+            positions[i] = x[i];
+        }
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Positions:" << std::endl;
+            for (std::vector<double>::iterator it = positions.begin();
+                 it != positions.end(); it += 3)
+            {
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+            }
+        }
+
+        // compute energy and forces using all MPI tasks
+        // expect positions to be replicated on all MPI tasks
+        std::vector<double> forces(n_);
+        double fx
+            = mgmol_->evaluateEnergyAndForces(positions, anumbers_, forces);
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Energy: " << fx << std::endl;
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                double norm = 0.;
+                for (int i = 0; i < 3; i++)
+                {
+                    double val = *(it + i);
+                    std::cout << "    " << val;
+                    norm += val * val;
+                }
+                std::cout << "  norm: " << std::sqrt(norm);
+                std::cout << std::endl;
+            }
+        }
+
+        // set gradient to negative forces
+        for (int i = 0; i < n_; i++)
+        {
+            grad[i] = -1. * forces[i];
+        }
+        return fx;
+    }
+};
+
+int main(int argc, char** argv)
+{
+    int mpirc = MPI_Init(&argc, &argv);
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Initialization failed!!!" << std::endl;
+        MPI_Abort(MPI_COMM_WORLD, 0);
+    }
+
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    /*
+     * Initialize general things, like magma, openmp, IO, ...
+     */
+    mgmol_init(comm);
+
+    /*
+     * read runtime parameters
+     */
+    std::string input_filename("");
+    std::string lrs_filename;
+    std::string constraints_filename("");
+
+    float total_spin = 0.;
+    bool with_spin   = false;
+
+    po::variables_map vm;
+
+    // read from PE0 only
+    if (MPIdata::onpe0)
+    {
+        read_config(argc, argv, vm, input_filename, lrs_filename,
+            constraints_filename, total_spin, with_spin);
+    }
+
+    MGmol_MPI::setup(comm, std::cout, with_spin);
+    MGmol_MPI& mmpi      = *(MGmol_MPI::instance());
+    MPI_Comm global_comm = mmpi.commGlobal();
+
+    /*
+     * Setup control struct with run time parameters
+     */
+    Control::setup(global_comm, with_spin, total_spin);
+    Control& ct = *(Control::instance());
+
+    ct.setOptions(vm);
+
+    int ret = ct.checkOptions();
+    if (ret < 0) return ret;
+
+    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(lrs_filename);
+
+    // Enter main scope
+    {
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "Construct MGmol object..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+
+        MGmolInterface* mgmol;
+        if (ct.isLocMode())
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+        else
+            mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "MGmol setup..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+        mgmol->setup();
+
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "Setup done..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+
+        // here we just use the atomic positions read in and used
+        // to initialize MGmol
+        std::vector<double> positions;
+        mgmol->getAtomicPositions(positions);
+        std::vector<short> anumbers;
+        mgmol->getAtomicNumbers(anumbers);
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Positions:" << std::endl;
+            std::vector<short>::iterator ita = anumbers.begin();
+            for (std::vector<double>::iterator it = positions.begin();
+                 it != positions.end(); it += 3)
+            {
+                std::cout << *ita;
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+                ita++;
+            }
+        }
+
+        // Set up parameters
+        LBFGSpp::LBFGSParam<double> param;
+        param.epsilon        = 4e-4;
+        param.max_iterations = 100;
+
+        // Create solver and function object
+        LBFGSpp::LBFGSSolver<double> solver(param);
+        const int n = positions.size();
+        if (MPIdata::onpe0) std::cout << "n = " << n << std::endl;
+        MGmolEnergyAndForces fun(mgmol, n, anumbers);
+
+        // initial guess
+        VectorXd x = VectorXd::Zero(n);
+        int i      = 0;
+        for (auto& pos : positions)
+        {
+            x[i++] = pos;
+        }
+
+        double eks;
+        int niter = solver.minimize(fun, x, eks);
+
+        std::cout << niter << " iterations" << std::endl;
+
+        delete mgmol;
+
+    } // close main scope
+
+    mgmol_finalize();
+
+    mpirc = MPI_Finalize();
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
+    }
+
+    time_t tt;
+    time(&tt);
+    if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl;
+
+    return 0;
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 3e1188fa..0f1edb08 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -68,7 +68,7 @@ set(SOURCES
  GrassmanLineMinimization.cc 
  GrassmanCG.cc 
  GrassmanCGSparse.cc 
- LBFGS.cc 
+ MGmolLBFGS.cc
  IonicStepper.cc 
  Energy.cc 
  GramMatrix.cc 
diff --git a/src/MGmol.cc b/src/MGmol.cc
index d5daabe6..f400d6ed 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -37,12 +37,12 @@
 #include "Hamiltonian.h"
 #include "Ions.h"
 #include "KBPsiMatrixSparse.h"
-#include "LBFGS.h"
 #include "LocGridOrbitals.h"
 #include "LocalizationRegions.h"
 #include "MDfiles.h"
 #include "MGkernels.h"
 #include "MGmol.h"
+#include "MGmolLBFGS.h"
 #include "MLWFTransform.h"
 #include "MPIdata.h"
 #include "MasksSet.h"
diff --git a/src/LBFGS.cc b/src/MGmolLBFGS.cc
similarity index 86%
rename from src/LBFGS.cc
rename to src/MGmolLBFGS.cc
index 9154b7b3..3cbfedcc 100644
--- a/src/LBFGS.cc
+++ b/src/MGmolLBFGS.cc
@@ -7,7 +7,7 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-#include "LBFGS.h"
+#include "MGmolLBFGS.h"
 #include "Control.h"
 #include "Electrostatic.h"
 #include "MGmol.h"
@@ -16,7 +16,7 @@
 #include "ProjectedMatrices.h"
 
 template <class OrbitalsType>
-LBFGS<OrbitalsType>::LBFGS(OrbitalsType** orbitals, Ions& ions,
+MGmolLBFGS<OrbitalsType>::MGmolLBFGS(OrbitalsType** orbitals, Ions& ions,
     Rho<OrbitalsType>& rho, ConstraintSet& constraints,
     std::shared_ptr<LocalizationRegions> lrs, ClusterOrbitals* local_cluster,
     MasksSet& masks, MasksSet& corrmasks, Electrostatic& electrostat,
@@ -37,7 +37,7 @@ LBFGS<OrbitalsType>::LBFGS(OrbitalsType** orbitals, Ions& ions,
 }
 
 template <class OrbitalsType>
-void LBFGS<OrbitalsType>::setup(const double dt)
+void MGmolLBFGS<OrbitalsType>::setup(const double dt)
 {
     if (lrs_)
         ref_lrs_ = std::shared_ptr<LocalizationRegions>(
@@ -77,14 +77,14 @@ void LBFGS<OrbitalsType>::setup(const double dt)
 }
 
 template <class OrbitalsType>
-LBFGS<OrbitalsType>::~LBFGS()
+MGmolLBFGS<OrbitalsType>::~MGmolLBFGS()
 {
     delete vh_init_;
     delete stepper_;
 }
 
 template <class OrbitalsType>
-void LBFGS<OrbitalsType>::reset(const double dt)
+void MGmolLBFGS<OrbitalsType>::reset(const double dt)
 {
     delete vh_init_;
     delete stepper_;
@@ -93,7 +93,7 @@ void LBFGS<OrbitalsType>::reset(const double dt)
 }
 
 template <class OrbitalsType>
-int LBFGS<OrbitalsType>::quenchElectrons(const int itmax, double& etot)
+int MGmolLBFGS<OrbitalsType>::quenchElectrons(const int itmax, double& etot)
 {
     etot_i_[0] = etot_i_[1];
     etot_i_[1] = etot_i_[2];
@@ -104,7 +104,7 @@ int LBFGS<OrbitalsType>::quenchElectrons(const int itmax, double& etot)
 }
 
 template <class OrbitalsType>
-void LBFGS<OrbitalsType>::updateRefs()
+void MGmolLBFGS<OrbitalsType>::updateRefs()
 {
     if (!stepper_->check_last_step_accepted())
     {
@@ -128,7 +128,7 @@ void LBFGS<OrbitalsType>::updateRefs()
 }
 
 template <>
-void LBFGS<LocGridOrbitals>::updateRefMasks()
+void MGmolLBFGS<LocGridOrbitals>::updateRefMasks()
 {
     Control& ct = *(Control::instance());
 
@@ -147,12 +147,12 @@ void LBFGS<LocGridOrbitals>::updateRefMasks()
 }
 
 template <>
-void LBFGS<ExtendedGridOrbitals>::updateRefMasks()
+void MGmolLBFGS<ExtendedGridOrbitals>::updateRefMasks()
 {
 }
 
 template <class OrbitalsType>
-void LBFGS<OrbitalsType>::setQuenchTol() const
+void MGmolLBFGS<OrbitalsType>::setQuenchTol() const
 {
     Control& ct = *(Control::instance());
     ct.conv_tol = 0.1 * stepper_->etol();
@@ -163,7 +163,7 @@ void LBFGS<OrbitalsType>::setQuenchTol() const
 }
 
 template <class OrbitalsType>
-void LBFGS<OrbitalsType>::updatePotAndMasks()
+void MGmolLBFGS<OrbitalsType>::updatePotAndMasks()
 {
     if (!stepper_->check_last_step_accepted())
         electrostat_.setupInitialVh(*vh_init_);
@@ -172,10 +172,10 @@ void LBFGS<OrbitalsType>::updatePotAndMasks()
 }
 
 template <class OrbitalsType>
-bool LBFGS<OrbitalsType>::lbfgsLastStepNotAccepted() const
+bool MGmolLBFGS<OrbitalsType>::lbfgsLastStepNotAccepted() const
 {
     return !stepper_->check_last_step_accepted();
 }
 
-template class LBFGS<LocGridOrbitals>;
-template class LBFGS<ExtendedGridOrbitals>;
+template class MGmolLBFGS<LocGridOrbitals>;
+template class MGmolLBFGS<ExtendedGridOrbitals>;
diff --git a/src/LBFGS.h b/src/MGmolLBFGS.h
similarity index 92%
rename from src/LBFGS.h
rename to src/MGmolLBFGS.h
index a085af05..c8348d24 100644
--- a/src/LBFGS.h
+++ b/src/MGmolLBFGS.h
@@ -27,7 +27,7 @@ class MGmol;
 class KBPsiMatrixInterface;
 
 template <class OrbitalsType>
-class LBFGS : public IonicAlgorithm<OrbitalsType>
+class MGmolLBFGS : public IonicAlgorithm<OrbitalsType>
 {
 private:
     OrbitalsType** orbitals_;
@@ -54,12 +54,12 @@ class LBFGS : public IonicAlgorithm<OrbitalsType>
     void setup(const double dt);
 
 public:
-    LBFGS(OrbitalsType** orbitals, Ions& ions, Rho<OrbitalsType>& rho,
+    MGmolLBFGS(OrbitalsType** orbitals, Ions& ions, Rho<OrbitalsType>& rho,
         ConstraintSet& constraints, std::shared_ptr<LocalizationRegions> lrs,
         ClusterOrbitals* local_cluster, MasksSet& masks, MasksSet& corrmasks,
         Electrostatic& electrostat, const double dt, MGmol<OrbitalsType>&);
 
-    ~LBFGS() override;
+    ~MGmolLBFGS() override;
 
     void reset(const double dt);
 
diff --git a/src/MGmol_NEB.cc b/src/MGmol_NEB.cc
index 5008ea30..b0f70e6f 100644
--- a/src/MGmol_NEB.cc
+++ b/src/MGmol_NEB.cc
@@ -9,8 +9,8 @@
 
 #include "DFTsolver.h"
 #include "FIRE.h"
-#include "LBFGS.h"
 #include "MGmol.h"
+#include "MGmolLBFGS.h"
 #include "ProjectedMatricesInterface.h"
 
 using namespace std;
@@ -34,8 +34,8 @@ void MGmol<T>::geomOptimSetup()
     switch (ct.AtomsDynamic())
     {
         case AtomsDynamicType::LBFGS:
-            geom_optimizer_ = new LBFGS<T>(&current_orbitals_, *ions_, *rho_,
-                *constraints_, *lrs_, local_cluster_, *currentMasks_,
+            geom_optimizer_ = new MGmolLBFGS<T>(&current_orbitals_, *ions_,
+                *rho_, *constraints_, *lrs_, local_cluster_, *currentMasks_,
                 *corrMasks_, *electrostat_, ct.dt, *this);
             break;
 
diff --git a/src/lbfgsrlx.cc b/src/lbfgsrlx.cc
index 6698eef0..b66305a1 100644
--- a/src/lbfgsrlx.cc
+++ b/src/lbfgsrlx.cc
@@ -13,10 +13,10 @@
 #include "Electrostatic.h"
 #include "Energy.h"
 #include "Ions.h"
-#include "LBFGS.h"
 #include "LBFGS_IonicStepper.h"
 #include "LocalizationRegions.h"
 #include "MGmol.h"
+#include "MGmolLBFGS.h"
 #include "MGmol_blas1.h"
 #include "MPIdata.h"
 #include "MasksSet.h"
@@ -35,7 +35,7 @@ void MGmol<OrbitalsType>::lbfgsrlx(OrbitalsType** orbitals, Ions& ions)
 {
     Control& ct = *(Control::instance());
 
-    LBFGS<OrbitalsType> lbfgs(orbitals, ions, *rho_, *constraints_, lrs_,
+    MGmolLBFGS<OrbitalsType> lbfgs(orbitals, ions, *rho_, *constraints_, lrs_,
         local_cluster_.get(), *currentMasks_, *corrMasks_, *electrostat_, ct.dt,
         *this);