Add new functionality to compute energy and forces (#267)

* use specified initial conditions for wavefunctions

src/Control.cc

@@ -997,42 +997,6 @@ void Control::readRestartInfo(std::ifstream* tfile)
     }
 }
 
-void Control::readRestartOutputInfo(std::ifstream* tfile)
-{
-    const std::string zero = "0";
-    const std::string one  = "1";
-    if (tfile != nullptr)
-    {
-        // Read in the output restart filename
-        std::string filename;
-        (*tfile) >> filename;
-        // int dpcs_chkpoint=0;
-        if (zero.compare(filename) == 0) // no restart dump
-            out_restart_info = 0;
-        else
-        {
-            if (one.compare(filename) == 0)
-            { // automatic naming of dump
-                filename                         = "snapshot";
-                out_restart_file_naming_strategy = 1;
-            }
-            (*tfile) >> out_restart_info;
-            (*tfile) >> out_restart_file_type;
-            //(*tfile)>>dpcs_chkpoint;
-            // timeout_.set(dpcs_chkpoint);
-        }
-
-        out_restart_file.assign(run_directory_);
-        out_restart_file.append("/");
-        out_restart_file.append(filename);
-        (*MPIdata::sout) << "Output restart file: " << out_restart_file
-                         << " with info level " << out_restart_info
-                         << std::endl;
-        //(*MPIdata::sout)<<"Time for DPCS checkpoint:
-        //"<<dpcs_chkpoint<<"[s]"<<endl;
-    }
-}
-
 int Control::setPreconditionerParameters(const short type, const float factor,
     const bool project_out, const short nlevels, const float fgrid_hmax)
 {

src/Control.h

@@ -284,7 +284,6 @@ class Control
     void adjust();
     int checkState();
     void readRestartInfo(std::ifstream* tfile);
-    void readRestartOutputInfo(std::ifstream* tfile);
     int readThermostatInfo(std::ifstream* tfile);
     void printThermostatInfo(std::ostream& os) const;
     void setNumst(const short myspin, const int nval);

src/IonicAlgorithm.cc

@@ -98,7 +98,7 @@ void IonicAlgorithm<OrbitalsType>::init(HDFrestart* h5f_file)
 template <class OrbitalsType>
 int IonicAlgorithm<OrbitalsType>::quenchElectrons(const int itmax, double& etot)
 {
-    int ret = mgmol_strategy_.quench(*orbitals_, ions_, itmax, 0, etot);
+    int ret = mgmol_strategy_.quench(**orbitals_, ions_, itmax, 0, etot);
 
     return ret;
 }

src/LBFGS.cc

@@ -97,7 +97,7 @@ int LBFGS<OrbitalsType>::quenchElectrons(const int itmax, double& etot)
 {
     etot_i_[0] = etot_i_[1];
     etot_i_[1] = etot_i_[2];
-    int ret = mgmol_strategy_.quench(*orbitals_, ions_, itmax, 0, etot_i_[2]);
+    int ret = mgmol_strategy_.quench(**orbitals_, ions_, itmax, 0, etot_i_[2]);
 
     etot = etot_i_[2];
     return ret;

src/MGmol.cc

@@ -521,7 +521,8 @@ void MGmol<OrbitalsType>::run()
     switch (ct.AtomsDynamic())
     {
         case AtomsDynamicType::Quench:
-            quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks);
+            quench(
+                *current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks);
 
             // Forces for the last states
             force(*current_orbitals_, *ions_);
@@ -1097,25 +1098,21 @@ void MGmol<OrbitalsType>::setup()
 }
 
 template <class OrbitalsType>
-void MGmol<OrbitalsType>::cleanup()
+void MGmol<OrbitalsType>::dumpRestart()
 {
-    closing_tm_.start();
-
-    Mesh* mymesh             = Mesh::instance();
-    const pb::PEenv& myPEenv = mymesh->peenv();
-    Control& ct              = *(Control::instance());
-
-    printTimers();
+    Control& ct = *(Control::instance());
 
     // Save data to restart file
-    if (ct.out_restart_info > 0 && !ct.AtomsMove())
+    if (ct.out_restart_info > 0)
     {
+        Mesh* mymesh           = Mesh::instance();
         const pb::Grid& mygrid = mymesh->grid();
         unsigned gdim[3] = { mygrid.gdim(0), mygrid.gdim(1), mygrid.gdim(2) };
+        const pb::PEenv& myPEenv = mymesh->peenv();
 
         // create restart file
         std::string filename(std::string(ct.out_restart_file));
-        filename += "0";
+        if (ct.out_restart_file_naming_strategy) filename += "0";
         HDFrestart h5restartfile(
             filename, myPEenv, gdim, ct.out_restart_file_type);
 
@@ -1125,6 +1122,22 @@ void MGmol<OrbitalsType>::cleanup()
         if (ierr < 0)
             os_ << "WARNING: writing restart data failed!!!" << std::endl;
     }
+}
+
+template <class OrbitalsType>
+void MGmol<OrbitalsType>::cleanup()
+{
+    closing_tm_.start();
+
+    Control& ct = *(Control::instance());
+
+    printTimers();
+
+    // Save data to restart file
+    if (!ct.AtomsMove())
+    {
+        dumpRestart();
+    }
 
     MPI_Barrier(comm_);
     closing_tm_.stop();
@@ -1421,6 +1434,14 @@ template <class OrbitalsType>
 double MGmol<OrbitalsType>::evaluateEnergyAndForces(
     const std::vector<double>& tau, std::vector<short>& atnumbers,
     std::vector<double>& forces)
+{
+    return evaluateEnergyAndForces(current_orbitals_, tau, atnumbers, forces);
+}
+
+template <class OrbitalsType>
+double MGmol<OrbitalsType>::evaluateEnergyAndForces(Orbitals* orbitals,
+    const std::vector<double>& tau, std::vector<short>& atnumbers,
+    std::vector<double>& forces)
 {
     assert(tau.size() == 3 * atnumbers.size());
 
@@ -1430,10 +1451,11 @@ double MGmol<OrbitalsType>::evaluateEnergyAndForces(
 
     moveVnuc(*ions_);
 
-    double eks = 0.;
-    quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks);
+    double eks              = 0.;
+    OrbitalsType* dorbitals = dynamic_cast<OrbitalsType*>(orbitals);
+    quench(*dorbitals, *ions_, ct.max_electronic_steps, 20, eks);
 
-    force(*current_orbitals_, *ions_);
+    force(*dorbitals, *ions_);
 
     ions_->getForces(forces);
 

src/MGmol.h

@@ -42,7 +42,6 @@ class IonicAlgorithm;
 #include "AOMMprojector.h"
 #include "ClusterOrbitals.h"
 #include "DMStrategy.h"
-#include "Energy.h"
 #include "ExtendedGridOrbitals.h"
 #include "Forces.h"
 #include "Ions.h"
@@ -131,7 +130,7 @@ class MGmol : public MGmolInterface
         KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix<DISTMATDTYPE>& hij);
     void computeHnlPhiAndAdd2HPhi(Ions& ions, OrbitalsType& phi,
         OrbitalsType& hphi, const KBPsiMatrixSparse* const kbpsi);
-    int dumprestartFile(OrbitalsType** orbitals, Ions& ions,
+    int dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
         Rho<OrbitalsType>& rho, const bool write_extrapolated_wf,
         const short count);
 
@@ -190,6 +189,9 @@ class MGmol : public MGmolInterface
     double evaluateEnergyAndForces(const std::vector<double>& tau,
         std::vector<short>& atnumbers, std::vector<double>& forces);
 
+    double evaluateEnergyAndForces(Orbitals*, const std::vector<double>& tau,
+        std::vector<short>& atnumbers, std::vector<double>& forces);
+
     /*
      * get internal atomic positions
      */
@@ -247,7 +249,7 @@ class MGmol : public MGmolInterface
 
     void update_pot(const pb::GridFunc<POTDTYPE>& vh_init, const Ions& ions);
     void update_pot(const Ions& ions);
-    int quench(OrbitalsType* orbitals, Ions& ions, const int max_steps,
+    int quench(OrbitalsType& orbitals, Ions& ions, const int max_steps,
         const int iprint, double& last_eks);
     int outerSolve(OrbitalsType& orbitals, OrbitalsType& work_orbitals,
         Ions& ions, const int max_steps, const int iprint, double& last_eks);
@@ -320,7 +322,17 @@ class MGmol : public MGmolInterface
         forces_->force(orbitals, ions);
     }
 
+    /*
+     * simply dump current state
+     */
+    void dumpRestart();
+
     void loadRestartFile(const std::string filename);
+
+    std::shared_ptr<ProjectedMatricesInterface> getProjectedMatrices()
+    {
+        return proj_matrices_;
+    }
 };
 // Instantiate static variables here to avoid clang warnings
 template <class OrbitalsType>

src/MGmolInterface.h

@@ -28,8 +28,15 @@ class MGmolInterface
     virtual double evaluateEnergyAndForces(const std::vector<double>& tau,
         std::vector<short>& atnumbers, std::vector<double>& forces)
         = 0;
+    virtual double evaluateEnergyAndForces(Orbitals*,
+        const std::vector<double>& tau, std::vector<short>& atnumbers,
+        std::vector<double>& forces)
+        = 0;
     virtual void getAtomicPositions(std::vector<double>& tau) = 0;
     virtual void getAtomicNumbers(std::vector<short>& an)     = 0;
+    virtual std::shared_ptr<ProjectedMatricesInterface> getProjectedMatrices()
+        = 0;
+    virtual void dumpRestart() = 0;
 };
 
 #endif

src/Orbitals.h

@@ -26,6 +26,8 @@ class Orbitals
 
     Orbitals() { iterative_index_ = -10; }
 
+    virtual ~Orbitals(){};
+
     Orbitals(const Orbitals& A, const bool copy_data)
     {
         if (copy_data)

src/md.cc

@@ -226,7 +226,7 @@ void checkMaxForces(const std::vector<double>& fion,
 }
 
 template <class OrbitalsType>
-int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
+int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
     Rho<OrbitalsType>& rho, const bool write_extrapolated_wf, const short count)
 {
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
@@ -255,7 +255,7 @@ int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
     {
         if (onpe0)
             (*MPIdata::serr)
-                << "dumprestartFile: cannot write ...previous_orbitals..."
+                << "dumpMDrestartFile: cannot write ...previous_orbitals..."
                 << std::endl;
         return ierr;
     }
@@ -269,7 +269,7 @@ int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
         if (ierr < 0)
         {
             if (onpe0)
-                (*MPIdata::serr) << "dumprestartFile: cannot write "
+                (*MPIdata::serr) << "dumpMDrestartFile: cannot write "
                                     "...ExtrapolatedFunction..."
                                  << std::endl;
             return ierr;
@@ -282,7 +282,7 @@ int MGmol<OrbitalsType>::dumprestartFile(OrbitalsType** orbitals, Ions& ions,
     {
         if (onpe0)
             (*MPIdata::serr)
-                << "dumprestartFile: cannot close file..." << std::endl;
+                << "dumpMDrestartFile: cannot close file..." << std::endl;
         return ierr;
     }
 
@@ -398,7 +398,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
     if (ct.restart_info < 3)
     {
         double eks = 0.;
-        quench(*orbitals, ions, ct.max_electronic_steps, 20, eks);
+        quench(**orbitals, ions, ct.max_electronic_steps, 20, eks);
     }
 
     ct.max_changes_pot = 0;
@@ -426,7 +426,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         bool last_move_is_small = true;
         do
         {
-            retval = quench(*orbitals, ions, ct.max_electronic_steps, 0, eks);
+            retval = quench(**orbitals, ions, ct.max_electronic_steps, 0, eks);
 
             // update localization regions
             if (ct.adaptiveLRs())
@@ -614,12 +614,12 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
                     while (ierr < 0 && count < DUMP_MAX_NUM_TRY)
                     {
                         dump_tm_.start();
-                        ierr = dumprestartFile(
+                        ierr = dumpMDrestartFile(
                             orbitals, ions, *rho_, extrapolated_flag, count);
                         dump_tm_.stop();
                         if (onpe0 && ierr < 0 && count < (DUMP_MAX_NUM_TRY - 1))
                             std::cout
-                                << "dumprestartFile() failed... try again..."
+                                << "dumpMDrestartFile() failed... try again..."
                                 << std::endl;
                         if (ierr < 0) sleep(1.);
                         count++;
@@ -640,12 +640,12 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         while (ierr < 0 && count < DUMP_MAX_NUM_TRY)
         {
             dump_tm_.start();
-            ierr = dumprestartFile(
+            ierr = dumpMDrestartFile(
                 orbitals, ions, *rho_, extrapolated_flag, count);
             dump_tm_.stop();
 
             if (onpe0 && ierr < 0 && count < (DUMP_MAX_NUM_TRY - 1))
-                std::cout << "dumprestartFile() failed... try again..."
+                std::cout << "dumpMDrestartFile() failed... try again..."
                           << std::endl;
             if (ierr < 0) sleep(1.);
             count++;