restart afqmc calcualtion

sanshar · Dec 18, 2023 · 35e196b · 35e196b
1 parent 4fcc287
commit 35e196b
Show file tree

Hide file tree

Showing 3 changed files with 103 additions and 2 deletions.
diff --git a/DQMC/MixedEstimator.cpp b/DQMC/MixedEstimator.cpp
@@ -5,6 +5,7 @@
 #include "DQMCStatistics.h"
 #include "MixedEstimator.h"
 #include <boost/format.hpp>
+#include <highfive/H5Easy.hpp>
 
 using namespace std;
 using namespace Eigen;
@@ -489,7 +490,86 @@ void calcMixedEstimatorLongProp(Wavefunction& waveLeft, Wavefunction& waveRight,
   double weightCap = 0.;
   if (schd.weightCap > 0) weightCap = schd.weightCap;
   else weightCap = std::max(100., walkers.size() / 10.);
-  for (int step = 1; step < nsweeps * nsteps; step++) {
+
+  // read restart information
+  if (schd.dqmcRestart) {
+    int nchol;
+    int read_sweeps;
+    vector<complex<double>> saved_walkers;
+    vector<complex<double>> saved_overlaps;
+    vector<double> saved_weights;
+    vector<double> saved_energies;
+    double totalCumulativeWeight;
+    if (commrank == 0) {
+      H5Easy::File chkfile("afqmc.h5", H5Easy::File::ReadOnly);
+      measureCounter = H5Easy::load<int>(chkfile, "block");
+      nchol = H5Easy::load<int>(chkfile, "nchol");
+      saved_walkers = H5Easy::load<std::vector<complex<double>>>(chkfile, "walkers");
+      saved_overlaps = H5Easy::load<std::vector<complex<double>>>(chkfile, "overlaps");
+      saved_weights = H5Easy::load<std::vector<double>>(chkfile, "weights");
+      saved_energies = H5Easy::load<std::vector<double>>(chkfile, "energies");
+      totalCumulativeWeight = H5Easy::load<double>(chkfile, "totalWeights");
+      eshift = H5Easy::load<double>(chkfile, "eshift");
+      eEstimate = H5Easy::load<double>(chkfile, "eEstimate");
+      read_sweeps = saved_weights.size();
+      cout << saved_walkers.size() << " " << saved_overlaps.size() << " " << saved_weights.size() << endl;
+    }
+    MPI_Barrier(MPI_COMM_WORLD);
+    int matSize;
+    if (schd.soc) matSize = 2 * ham.norbs * ham.nelec;
+    else if (walker.szQ) matSize = ham.norbs * ham.nelec;
+    else matSize = ham.norbs * (ham.nalpha  + ham.nbeta);
+    vector<complex<double>> serial(nwalk * matSize, complex<double>(0., 0.)), serialw(matSize, complex<double>(0., 0.)), overlaps(nwalk, complex<double>(0., 0.));
+    MPI_Bcast(&read_sweeps, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Bcast(&eshift, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+    MPI_Bcast(&eEstimate, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+    MPI_Bcast(&totalCumulativeWeight, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+    MPI_Bcast(&measureCounter, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    saved_walkers.resize(nwalk*matSize);
+    saved_overlaps.resize(nwalk);
+    saved_weights.resize(read_sweeps);
+    saved_energies.resize(read_sweeps);
+    MPI_Bcast(saved_energies.data(),
+              read_sweeps,
+              MPI_DOUBLE,
+              0,
+              MPI_COMM_WORLD);
+    MPI_Bcast(saved_walkers.data(), nwalk * matSize, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+    MPI_Bcast(saved_overlaps.data(), nwalk, MPI_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+    MPI_Bcast(saved_weights.data(), read_sweeps, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+    // deserialize
+    for (int w = 0; w < walkers.size(); w++) {
+      for (int i = 0; i < matSize; i++) {serialw[i] = saved_walkers[w * matSize + i];
+      cout << serialw[i] << endl;}
+      walkers[w].setDet(serialw, overlaps[w]);
+      weights[w] = totalCumulativeWeight / commsize / nwalk;
+      cout << walkers[w].trialOverlap << endl;
+      walkers[w].overlap(waveLeft);
+    }
+
+    for (int i = 0; i < read_sweeps; i++) {
+      totalEnergies(i) = saved_energies[i];
+      totalWeights(i) = saved_weights[i];
+      if (i < schd.burnIter) {
+        averageNumEql += totalEnergies(i)*totalWeights(i);
+        averageDenomEql += totalWeights(i);
+        averageEnergyEql = averageNumEql / averageDenomEql;
+        afqmcFile << boost::format(" %5d     %.3e       %.5s     %.5e      %.9e         %7s         %.2e \n") % i % 0.0 % '-' % totalWeights(i) % totalEnergies(i) % '-' % (getTime()-calcInitTime); 
+        afqmcFile.flush();
+      }
+      else {
+        averageNum += totalEnergies(i)*totalWeights(i);
+        averageDenom += totalWeights(i);
+        averageEnergy = averageNum / averageDenom;
+        afqmcFile << boost::format(" %5d     %.3e       %.5s     %.5e      %.9e        %.9e        %.2e \n") % i % 0.0 % '-' % totalWeights(i) % totalEnergies(i) % averageEnergy % (getTime() - calcInitTime); 
+        afqmcFile.flush();
+      }
+    }
+    totalWeight = totalCumulativeWeight / commsize;
+    MPI_Barrier(MPI_COMM_WORLD);
+  }
+  for (int step = 1+nsteps*(measureCounter-1); step < (nsweeps-1) * nsteps+1; step++) {
     // average before eql
     if (step * dt < 10.) averageEnergy = averageEnergyEql;
 
@@ -511,7 +591,6 @@ void calcMixedEstimatorLongProp(Wavefunction& waveLeft, Wavefunction& waveRight,
     MPI_Allreduce(MPI_IN_PLACE, &totalWeight, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
     //eshift = averageEnergy - 0.1 * log(totalWeight/(nwalk * commsize)) / dt;
     eshift = eEstimate - 0.1 * log(totalWeight/(nwalk * commsize)) / dt;
-
     // orthogonalize for stability
     if (step % orthoSteps == 0) {
       for (int w = 0; w < walkers.size(); w++) walkers[w].orthogonalize();
@@ -644,6 +723,22 @@ void calcMixedEstimatorLongProp(Wavefunction& waveLeft, Wavefunction& waveRight,
           for (int i = 0; i < matSize; i++) newSerialGather[w * matSize + i] = serialGather[index * matSize + i];
           newOverlapsGather[w] = overlapsGather[index];
         }
+        int block = step / nsteps;
+        if (block % 100 == 0) {
+          H5Easy::File chkfile("afqmc.h5", H5Easy::File::Overwrite);
+          if (schd.writeWalkers) {
+            H5Easy::dump(chkfile, "overlaps", newOverlapsGather, H5Easy::DumpMode::Overwrite);
+            H5Easy::dump(chkfile, "walkers", newSerialGather, H5Easy::DumpMode::Overwrite);
+            H5Easy::dump(chkfile, "totalWeights", totalCumulativeWeight, H5Easy::DumpMode::Overwrite);
+          }
+          H5Easy::dump(chkfile, "weights", totalWeights, H5Easy::DumpMode::Overwrite);
+          H5Easy::dump(chkfile, "energies", totalEnergies, H5Easy::DumpMode::Overwrite);
+          H5Easy::dump(chkfile, "nchol", nchol, H5Easy::DumpMode::Overwrite);
+          H5Easy::dump(chkfile, "block", block, H5Easy::DumpMode::Overwrite);
+          H5Easy::dump(chkfile, "averageEnergy", block, H5Easy::DumpMode::Overwrite);
+          H5Easy::dump(chkfile, "eEstimate", eEstimate, H5Easy::DumpMode::Overwrite);
+          H5Easy::dump(chkfile, "eshift", eshift, H5Easy::DumpMode::Overwrite);
+        }
       }
 
       // scatter

diff --git a/utils/input.cpp b/utils/input.cpp
@@ -243,6 +243,7 @@ void readInput(string inputFile, schedule& schd, bool print) {
     else if (method == "lm") schd.method = linearmethod;
     schd.restart = input.get("optimizer.restart", false);
     schd.fullRestart = input.get("optimizer.fullRestart", false);
+    schd.dqmcRestart = input.get("sampling.restart", false);
     child = input.get_child_optional("sampling.maxIter"); //to ensure maxiter is not reassigned
     if (!child) schd.maxIter = input.get("optimizer.maxIter", 50);
     schd.avgIter = input.get("optimizer.avgIter", 0);
@@ -279,6 +280,7 @@ void readInput(string inputFile, schedule& schd, bool print) {
     // dqmc
     schd.writeOneRDM = input.get("print.writeOneRDM", false); // also used in vmc
     schd.writeTwoRDM = input.get("print.writeTwoRDM", false); // used in vmc
+    schd.writeWalkers = input.get("print.writeWalkers", false);
     schd.scratchDir = input.get("print.scratchDir", "./");
 
 

diff --git a/utils/input.h b/utils/input.h
@@ -45,6 +45,7 @@ struct schedule {
     ar & restart & deterministic
       & tol & correlatorFiles
       & fullRestart
+      & dqmcRestart
       & wavefunctionType
       & integralsFile
       & ghfDets
@@ -162,6 +163,7 @@ struct schedule {
       & soc
       & writeOneRDM
       & writeTwoRDM
+      & writeWalkers
       & intType
       & scratchDir
       & weightCap
@@ -207,6 +209,7 @@ struct schedule {
 //General options
   bool restart;                          //option to restart calculation
   bool fullRestart;                          //option to restart calculation
+  bool dqmcRestart;                          //option to restart dqmc calculation
   bool deterministic;                    //Performs a deterministic calculation   
   int printLevel;                        // How much stuff to print
   bool expCorrelator;                    // exponential correlator parameters, to enforce positivity
@@ -402,6 +405,7 @@ struct schedule {
   bool soc;
   bool writeOneRDM;
   bool writeTwoRDM;
+  bool writeWalkers;
   std::string intType;
   std::string scratchDir;
   double weightCap;