LoadOrbitalsFromRestartFile and miscellaneous changes (#226)

* Signal.h -> mgmol_Signal.h

* Carbyne example files.

* fixed FindSCALAPACK.cmake to search for default paths

* MGmol::loadOrbitalFromRestartFile

* nullptr initialization

cmake_modules/FindSCALAPACK.cmake

@@ -7,7 +7,7 @@ if(DEFINED ENV{SCALAPACK_ROOT})
 endif(DEFINED ENV{SCALAPACK_ROOT})
 
 if(SCALAPACK_ROOT)
-  set(_SCALAPACK_SEARCH_DIR ${SCALAPACK_ROOT})
+  set(_SCALAPACK_SEARCH_DIR ${SCALAPACK_ROOT} ${SCALAPACK_ROOT}/lib/intel64)
   list(APPEND _SCALAPACK_SEARCHES ${_SCALAPACK_SEARCH_DIR})
 endif(SCALAPACK_ROOT)
 
@@ -29,13 +29,28 @@ if(NOT SCALAPACK_LIBRARY)
   endforeach()
 endif()
 
-unset(SCALAPACK_NAMES)
-
 mark_as_advanced(SCALAPACK_LIBRARY SCALAPACK_INCLUDE_DIR)
 
 include(FindPackageHandleStandardArgs)
 FIND_PACKAGE_HANDLE_STANDARD_ARGS(SCALAPACK REQUIRED_VARS SCALAPACK_LIBRARY)
 
+# Search for some default library paths
+if (NOT SCALAPACK_FOUND)
+  find_library(SCALAPACK_LIBRARY
+    NAMES ${SCALAPACK_NAMES}
+    PATHS /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib
+    /opt/local/lib /opt/sw/lib /sw/lib
+    ENV LD_LIBRARY_PATH
+    ENV DYLD_FALLBACK_LIBRARY_PATH
+    ENV DYLD_LIBRARY_PATH
+    ENV SCALAPACKDIR
+    ENV BLACSDIR)
+
+  FIND_PACKAGE_HANDLE_STANDARD_ARGS(SCALAPACK REQUIRED_VARS SCALAPACK_LIBRARY)
+endif()
+
+unset(SCALAPACK_NAMES)
+
 if(SCALAPACK_FOUND)
   # Only Intel's scalapack requires an include directory
   if(SCALAPACK_INCLUDE_DIR)

examples/Carbyne/carbyne.cfg

@@ -0,0 +1,26 @@
+verbosity=2
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx= 96
+ny= 96
+nz= 192
+[Domain]
+ox= -10.
+oy= -10.
+oz= -20.
+lx= 20.
+ly= 20.
+lz= 40.
+[Potentials]
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+pseudopotential=pseudo.C_ONCV_PBE_SG15
+[Run]
+type=QUENCH
+[Quench]
+max_steps=200
+atol=1.e-8
+[Orbitals]
+initial_type=Fourier
+[Restart]
+output_level=3

examples/Carbyne/carbyne.in

@@ -0,0 +1,14 @@
+H00 1  -0.0000   -0.0000   15.2674
+C01 2  -0.0000    0.0000   13.2519
+C02 2  -0.0000    0.0000   10.9495
+C03 2  -0.0000   -0.0000    8.4221
+C04 2   0.0000    0.0000    6.0897
+C05 2  -0.0000    0.0000    3.5892
+C06 2  -0.0000   -0.0000    1.2470
+C07 2   0.0000   -0.0000   -1.2469
+C08 2   0.0000   -0.0000   -3.5891
+C09 2  -0.0000   -0.0000   -6.0897
+C10 2  -0.0000    0.0000   -8.4221
+C11 2   0.0000   -0.0000  -10.9495
+C12 2   0.0000    0.0000  -13.2520
+H13 1   0.0000    0.0000  -15.2675

src/MGmol.cc

@@ -106,7 +106,7 @@ extern Timer ions_setupInteractingIons_tm;
 extern Timer ions_setup_tm;
 extern Timer updateCenters_tm;
 
-#include "Signal.h"
+#include "mgmol_Signal.h"
 std::set<int> Signal::recv_;
 
 template <class OrbitalsType>

src/MGmol.h

@@ -102,7 +102,7 @@ class MGmol : public MGmolInterface
     double total_energy_;
     ConstraintSet* constraints_;
 
-    OrbitalsExtrapolation<OrbitalsType>* orbitals_extrapol_;
+    OrbitalsExtrapolation<OrbitalsType>* orbitals_extrapol_ = nullptr;
 
     float md_time_;
     int md_iteration_;
@@ -301,6 +301,8 @@ class MGmol : public MGmolInterface
     {
         forces_->force(orbitals, ions);
     }
+
+    OrbitalsType* loadOrbitalFromRestartFile(const std::string filename);
 };
 // Instantiate static variables here to avoid clang warnings
 template <class OrbitalsType>

src/OrbitalsExtrapolation.h

@@ -44,7 +44,7 @@ class OrbitalsExtrapolation
     virtual short getNumOrbitalExtrapolations() { return 0; }
 
 protected:
-    OrbitalsType* orbitals_minus1_;
+    OrbitalsType* orbitals_minus1_ = nullptr;
 };
 
 #endif

src/md.cc

@@ -31,7 +31,7 @@
 #include "ProjectedMatricesMehrstellen.h"
 #include "ProjectedMatricesSparse.h"
 #include "Rho.h"
-#include "Signal.h"
+#include "mgmol_Signal.h"
 #include "SpreadsAndCenters.h"
 #include "tools.h"
 
@@ -674,5 +674,89 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
     delete orbitals_extrapol_;
 }
 
+template <class OrbitalsType>
+OrbitalsType* MGmol<OrbitalsType>::loadOrbitalFromRestartFile(const std::string filename)
+{
+    MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
+    Control& ct              = *(Control::instance());
+    Mesh* mymesh             = Mesh::instance();
+    const pb::PEenv& myPEenv = mymesh->peenv();
+
+    /* For now, we only consider double-precision hdf5 I/O. */
+    assert(ct.restart_info > 3);
+    assert((ct.AtomsDynamic() == AtomsDynamicType::MD) || (ct.AtomsDynamic() == AtomsDynamicType::Quench));
+
+    HDFrestart h5file(filename, myPEenv, ct.out_restart_file_type);
+    int ierr;
+
+    OrbitalsType *restart_orbitals = nullptr;
+
+    /* This corresponds to MGmol<OrbitalsType>::initial */
+    {
+        // Copy from current orbital, instead of constructing brand-new one
+        restart_orbitals = new OrbitalsType("ForLoading", *current_orbitals_, false);
+
+        /* This corresponds to MGmol<OrbitalsType>::read_restart_data */
+        {
+            ierr = restart_orbitals->read_func_hdf5(h5file);
+        }   // read_restart_data
+    }   // initial()
+
+    /* This corresponds to MGmol<OrbitalsType>::md */
+    {
+        int flag_extrapolated_data = 0;
+        if (onpe0)
+        {
+            flag_extrapolated_data
+                = h5file.dset_exists("ExtrapolatedFunction0000");
+            if (flag_extrapolated_data == 0)
+                flag_extrapolated_data
+                    = h5file.dset_exists("ExtrapolatedFunction0");
+        }
+        MPI_Bcast(&flag_extrapolated_data, 1, MPI_INT, 0, comm_);
+
+        /*
+            If extrapolated function exists,
+            then function is set as previous orbitals,
+            while the extrapolated function is set as the current orbitals.
+            This is how the restart file is saved via dumprestartFile.
+        */
+        if (flag_extrapolated_data)
+        {
+            orbitals_extrapol_ = OrbitalsExtrapolationFactory<OrbitalsType>::create(
+                                                                ct.WFExtrapolation());
+
+            if (onpe0) os_ << "Create new orbitals_minus1..." << std::endl;
+
+            orbitals_extrapol_->setupPreviousOrbitals(&restart_orbitals,
+                proj_matrices_, lrs_, local_cluster_, currentMasks_,
+                corrMasks_, h5file);
+        }
+
+        /* main workflow delete h5f_file_ here, meaning the loading is over. */
+    }   // md()
+
+    ierr = h5file.close();
+    mmpi.allreduce(&ierr, 1, MPI_MIN);
+    if (ierr < 0)
+    {
+        if (onpe0)
+            (*MPIdata::serr)
+                << "loadRestartFile: cannot close file..." << std::endl;
+        return nullptr;
+    }
+
+    /*
+        In returning restart_orbitals,
+        we hope that the wavefunctions in restart_orbitals are all set.
+        At least the following functions should return proper data loaded from the file:
+
+            restart_orbitals->getLocNumpt()
+            restart_orbitals->chromatic_number()
+            restart_orbitals->getPsi(idx)   (for int idx)
+    */
+    return restart_orbitals;
+}
+
 template class MGmol<LocGridOrbitals>;
 template class MGmol<ExtendedGridOrbitals>;

src/tools/Timeout.h

@@ -17,7 +17,7 @@
 #include <stdlib.h>
 
 #include "MPIdata.h"
-#include "Signal.h"
+#include "mgmol_Signal.h"
 
 #if PCS
 #include <csignal>

src/tools/Signal.h → src/tools/mgmol_Signal.h

@@ -15,7 +15,7 @@
 // access the flag set, reset flags, or interrogate flags.
 // The Signal class can be used in an application by declaring
 //
-// #include "Signal.h"
+// #include "mgmol_Signal.h"
 // set<int> Signal::recv_;
 //
 // A signal can be registered using, e.g.