Update ROMFPMD branch from release (#255)

* Reenable testShortSighted test (#248)

* Make new driver to check input (#247)

* clean up/reorganize main.cc
* use shared_ptr in class MGmol

* Add possible periodic dimensions to xyz2in.py (#249)

* Add possible periodic dimensions to xyz2in.py

* Remove unused/untested option extrapolateH (#250)

* Exit with failure if density off by more than 2% (#251)

* Exit with failure if density off by more than 2%
* adapt SiH4 test to catch that
* fix bug in DFTsolver that was leading to wrong density

* Example driver (#252)

* add example driver, showing use of MGmol as a force/energy computational engine

* clean up related functions in class Ions

* loadOrbitalsFromRestartFile -> loadRestartFile

* update the rom main driver.

* update rom main driver again

* temporary fix before the next PR

---------

Co-authored-by: Jean-Luc Fattebert <[email protected]>

CMakeLists.txt

@@ -178,7 +178,7 @@ if(${MGMOL_WITH_CLANG_FORMAT})
   find_package(CLANG_FORMAT)
   if(${CLANG_FORMAT_FOUND})
     message(STATUS "Indent with clang-format")
-    file(GLOB_RECURSE FORMAT_SOURCES src/*.cc src/*.h tests/*.cc tests/*.h)
+    file(GLOB_RECURSE FORMAT_SOURCES src/*.cc src/*.h tests/*.cc tests/*.h drivers/*.cc)
     add_custom_target(format
       COMMAND ${CLANG_FORMAT_EXECUTABLE} -i -style=file ${FORMAT_SOURCES}
       DEPENDS ${FORMAT_SOURCES})
@@ -271,5 +271,7 @@ link_libraries(${LIBROM_LIB})
 # add subdirectories for source files, tests
 add_subdirectory(src)
 
+add_subdirectory(drivers)
+
 add_subdirectory(tests)
 

drivers/CMakeLists.txt

@@ -0,0 +1,11 @@
+include_directories( ${CMAKE_SOURCE_DIR}/src )
+
+add_executable(check_input check_input.cc)
+
+add_executable(example1 example1.cc)
+
+target_include_directories(check_input PRIVATE ${Boost_INCLUDE_DIRS})
+target_include_directories(example1 PRIVATE ${Boost_INCLUDE_DIRS})
+
+target_link_libraries(check_input mgmol_src)
+target_link_libraries(example1 mgmol_src)

drivers/check_input.cc

@@ -0,0 +1,93 @@
+// 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
+
+#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 <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+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;
+
+    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 options 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();
+
+    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
+    {
+        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);
+
+        *MPIdata::sout << " Input parameters OK\n";
+
+        delete mgmol;
+    } // close main scope
+
+    mgmol_finalize();
+
+    mpirc = MPI_Finalize();
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
+    }
+
+    return 0;
+}

drivers/example1.cc

@@ -0,0 +1,167 @@
+// 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
+
+#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;
+
+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++;
+            }
+        }
+
+        // compute energy and forces using all MPI tasks
+        // expect positions to be replicated on all MPI tasks
+        std::vector<double> forces;
+        mgmol->evaluateEnergyAndForces(positions, anumbers, forces);
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << 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;
+}

src/CMakeLists.txt

@@ -191,13 +191,15 @@ install(TARGETS mgmol_src DESTINATION lib)
 add_executable(mgmol-opt main.cc)
 target_include_directories (mgmol-opt PRIVATE ${Boost_INCLUDE_DIRS})
 
-target_link_libraries(mgmol-opt mgmol_src ${link_libs})
-target_link_libraries(mgmol-opt ${SCALAPACK_LIBRARIES})
-target_link_libraries(mgmol-opt ${HDF5_LIBRARIES})
-target_link_libraries(mgmol-opt ${HDF5_HL_LIBRARIES})
-target_link_libraries(mgmol-opt ${BLAS_LIBRARIES})
-target_link_libraries(mgmol-opt ${LAPACK_LIBRARIES})
-target_link_libraries(mgmol-opt ${Boost_LIBRARIES})
+target_link_libraries(mgmol_src ${link_libs})
+target_link_libraries(mgmol_src ${SCALAPACK_LIBRARIES})
+target_link_libraries(mgmol_src ${LAPACK_LIBRARIES})
+target_link_libraries(mgmol_src ${BLAS_LIBRARIES})
+target_link_libraries(mgmol_src ${HDF5_LIBRARIES})
+target_link_libraries(mgmol_src ${HDF5_HL_LIBRARIES})
+target_link_libraries(mgmol_src ${Boost_LIBRARIES})
+
+target_link_libraries(mgmol-opt mgmol_src)
 if (${OPENMP_CXX_FOUND})
   target_link_libraries(mgmol-opt OpenMP::OpenMP_CXX)
 endif()

src/DFTsolver.cc

@@ -339,13 +339,15 @@ int DFTsolver<OrbitalsType>::solve(OrbitalsType& orbitals,
             else
             {
                 bool updateDM = false;
+                bool updatedS = false;
                 if (!ct.fullyOccupied())
                 {
                     orbitals.computeGramAndInvS();
                     dm_strategy_->dressDM();
                     updateDM = true;
+                    updatedS = true;
                 }
-                orbitals.orthonormalizeLoewdin(true, nullptr, updateDM);
+                orbitals.orthonormalizeLoewdin(updatedS, nullptr, updateDM);
 
                 orbitals_stepper_->restartMixing();
             }

src/Electrostatic.cc

@@ -36,7 +36,7 @@ Timer Electrostatic::solve_tm_("Electrostatic::solve");
 
 Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
     const double screening_const)
-    : laptype_(lap_type)
+    : laptype_(lap_type), poisson_solver_(nullptr)
 {
     assert(bcPoisson[0] >= 0);
     assert(bcPoisson[1] >= 0);
@@ -166,6 +166,8 @@ Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
 
 Electrostatic::~Electrostatic()
 {
+    assert(poisson_solver_ != nullptr);
+
     delete poisson_solver_;
     if (grhod_ != nullptr) delete grhod_;
     if (grhoc_ != nullptr) delete grhoc_;