add total cpu timings across processors; add docs sections and fixes.

Author: Jun Yang

dmrg.C

@@ -74,6 +74,7 @@ void license() {
 
 namespace SpinAdapted{
   Timer globaltimer(false);
+  double tcpu,twall,ecpu,ewall;
   bool DEBUGWAIT = false;
   bool DEBUG_MEMORY = false;
   bool restartwarm = false;
@@ -398,8 +399,10 @@ int calldmrg(char* input, char* output)
   }
 
   cout.rdbuf(backup);
-  pout << setprecision(3) <<"\n\n\t\t\t BLOCK CPU  Time (seconds): " << globaltimer.totalcputime() << endl;
-  pout << setprecision(3) <<"\t\t\t BLOCK Wall Time (seconds): " << globaltimer.totalwalltime() << endl;
+
+  tcpu=globaltimer.totalcputime();twall=globaltimer.totalwalltime();
+  pout << setprecision(3) <<"\n\n\t\t\t BLOCK CPU  Time (seconds): " << tcpu << endl;
+  pout << setprecision(3) <<"\t\t\t BLOCK Wall Time (seconds): " << twall << endl;
 
   return 0;
 }

docs/source/examples.rst

@@ -672,4 +672,20 @@ Energy extrapolation:
 
    Starting from *M=500*, use the largest discarded weights and associated sweep energies in the last sweep iteration of each *M* to make linear regression (see the figure above). The extrapolated DMRG sweep energy is -75.728557 a.u.
 
+Further Reading
+===============
+
+Some practical questions are often asked such as,
+
+* what sort of molecules can the DMRG be practically applied to?
+* what sort of accuracies can be obtained and at what cost? What are the typical sizes of systems (e.g. number
+  of active orbitals) that can be treated with practical computational resources?
+* how do we reason about the accuracy of DMRG calculations for dirent molecules?
+* how is a DMRG calculation best specified (e.g. in terms of starting orbitals and their order)?
+
+We provide answers in the following paper from both theoretical reasoning and numerical calculation by applying the DMRG 
+to a representative set of molecules. 
+The calculations we describe therein are all run in a completely black-box fashion using the default settings of our ``Block`` code.
+
+* R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang and G. K.-L. Chan, J. Chem. Phys. 142, 034102 (2015).
 

docs/source/index.rst

@@ -11,7 +11,7 @@ Welcome to BLOCK's documentation!
 
 `BLOCK` implements the density matrix renormalization group (DMRG) algorithm for quantum chemistry.
 The DMRG is a variational wavefunction method. Compared to other quantum chemical methods,
-it executionciently describes strong, multi-reference correlation in a large number of active orbitals (occupancies far from 0 or 2).
+it efficiently describes strong, multi-reference correlation in a large number of active orbitals (occupancies far from 0 or 2).
 The method is also provably optimal for correlation with a one-dimensional topology, that is,
 where orbitals are arranged with a chain- or ring-like connectivity. 
 However, with the possible exception of small molecules, for correlation that is dynamic in character, 
@@ -20,7 +20,7 @@ such as coupled cluster theory or multireference configuration interaction.
 We recommend the use of the DMRG in problems requiring active spaces too large for
 standard complete active space (CAS) techniques. Thus, if you are interested in:
 
-* a CAS-like treatment of low-lying eigenstates in problems with 16-40 active orbitals,
+* a CAS-like treatment of low-lying eigenstates in real problems with more than 50 active orbitals,
 * or, one-dimensional orbital topologies with up to 100 active orbitals,
 * and, standard chemical accuracy (1 kcal/mol in energy differences),
 

docs/source/overview.rst

@@ -10,7 +10,7 @@ Features
 * one-, two-, three- and four-particle density matrices,
 * one- and two-particle transition density matrices between two states,
 * perturbation methods including NEVPT2 and MPSPT,
-* interfaces to the Molpro, ORCA, Q-Chem and Molcas program packages. 
+* DMRG-SCF and/or DMRG-NEVPT2 interfaces to the Molpro, ORCA, Q-Chem and Molcas program packages. 
 
 Calling `Block` as an external function
 =======================================
@@ -34,8 +34,9 @@ use of `Block`:
 * G. K.-L. Chan and M. Head-Gordon, J. Chem. Phys. 116, 4462 (2002),
 * G. K.-L. Chan, J. Chem. Phys. 120, 3172 (2004),
 * D. Ghosh, J. Hachmann, T. Yanai, and G. K.-L. Chan, J. Chem. Phys., 128, 144117 (2008),
-* S. Sharma and G. K-.L. Chan, J. Chem. Phys. 136, 124121 (2012).
+* S. Sharma and G. K-.L. Chan, J. Chem. Phys. 136, 124121 (2012),
+* R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang and G. K.-L. Chan, J. Chem. Phys. 142, 034102 (2015).
 
 In addition, a useful list of DMRG references relevant to quantum chemistry can be found
-in the article above by Sharma and Chan.
+in the article above by Sharma and Chan. 
 

global.h

@@ -73,6 +73,7 @@ extern bool NonabelianSym;
 extern std::vector<int> NPROP;
 extern int PROPBITLEN;
 extern double NUMERICAL_ZERO;
+extern double tcpu,twall,ecpu,ewall;
 //extern ifstream* coutbuf;
 }
 #endif

modules/generate_blocks/sweep.C

@@ -182,8 +182,9 @@ double SweepGenblock::do_one(SweepParams &sweepParams, const bool &warmUp, const
 
   ++sweepParams.set_sweep_iter();
 
-  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << sweeptimer.elapsedcputime() << endl;
-  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << sweeptimer.elapsedwalltime() << endl;
+  ecpu = sweeptimer.elapsedcputime(); ewall = sweeptimer.elapsedwalltime();
+  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << ecpu << endl;
+  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << ewall << endl;
 
   return finalEnergy[0];
 }
@@ -252,8 +253,9 @@ void SweepGenblock::do_one(SweepParams &sweepParams, const bool &forward, int st
 
   ++sweepParams.set_sweep_iter();
 
-  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << sweeptimer.elapsedcputime() << endl;
-  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << sweeptimer.elapsedwalltime() << endl;
+  ecpu = sweeptimer.elapsedcputime(); ewall = sweeptimer.elapsedwalltime();
+  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << ecpu << endl;
+  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << ewall << endl;
 
 }
 

modules/npdm/fourpdm_container.C

@@ -74,7 +74,8 @@ void Fourpdm_container::save_npdms(const int& i, const int& j)
 #ifndef SERIAL
   world.barrier();
 #endif
-  p3out << "4PDM save full array time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+  ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+  p3out << "4PDM save full array time " << ewall << " " << ecpu << endl;
 
 }
 
@@ -283,18 +284,21 @@ void Fourpdm_container::save_spatial_npdm_binary(const int &i, const int &j)
       Sortpdm::partition_data<Sortpdm::index_element>((long)pow(dmrginp.last_site(),8),file,tmpfile);
       Sortpdm::externalsort<Sortpdm::index_element>(tmpfile,sortedfile,(long)pow(dmrginp.last_site(),8));
       world.barrier();
-      p3out << "4PDM parallel external sort time " << timer1.elapsedwalltime() << " " << timer1.elapsedcputime() << endl;
+      ecpu = timer1.elapsedcputime(); ewall = timer1.elapsedwalltime();
+      p3out << "4PDM parallel external sort time " << ewall << " " << ecpu << endl;
       Timer timer;
       Sortpdm::mergefile(sortedfile);
       world.barrier();
       if(mpigetrank()==0) boost::filesystem::rename(sortedfile,finalfile);
       boost::filesystem::remove(tmpfile);
-      p3out << "4PDM merge sorted file time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+      ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+      p3out << "4PDM merge sorted file time " << ewall << " " << ecpu << endl;
 #else
       Timer timer2;
       Sortpdm::externalsort<Sortpdm::index_element> (file,finalfile,(long)pow(dmrginp.last_site(),8));
       boost::filesystem::remove(file);
-      p3out << "4PDM external sort time " << timer2.elapsedwalltime() << " " << timer2.elapsedcputime() << endl;
+      ecpu = timer2.elapsedcputime();ewall=timer2.elapsedwalltime();
+      p3out << "4PDM external sort time " << ewall << " " << ecpu << endl;
 #endif
     }
   }

modules/npdm/npdm.C

@@ -195,7 +195,8 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP
   }
 
   //newSystem.transform_operators(rotateMatrix,rotateMatrixB);
-  p3out << "NPDM block and decimate and compute elements " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+  ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+  p3out << "NPDM block and decimate and compute elements " << ewall << " " << ecpu << endl;
 
 }
 
@@ -287,7 +288,8 @@ double npdm_do_one_sweep(Npdm_driver_base &npdm_driver, SweepParams &sweepParams
     ++sweepParams.set_block_iter();
     //sweepParams.savestate(forward, system.get_sites().size());
 
-    p3out << "NPDM do one site time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+    ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+    p3out << "NPDM do one site time " << ewall << " " << ecpu << endl;
   }
 
   //for(int j=0;j<nroots;++j)
@@ -308,8 +310,9 @@ double npdm_do_one_sweep(Npdm_driver_base &npdm_driver, SweepParams &sweepParams
   // Update the static number of iterations
   ++sweepParams.set_sweep_iter();
 
-  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << sweeptimer.elapsedcputime() << endl;
-  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << sweeptimer.elapsedwalltime() << endl;
+  ecpu = sweeptimer.elapsedcputime();ewall=sweeptimer.elapsedwalltime();
+  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << ecpu << endl;
+  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << ewall << endl;
 
   return finalEnergy[0];
 
@@ -400,7 +403,8 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
       SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[1]); //this will generate the cd operators                               
     else abort();
 		dmrginp.npdm_generate() = false;
-    p3out << "\t\t\t NPDM SweepGenblock time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+    ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+    p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl;
     dmrginp.set_fullrestart() = false;
 
     sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy;
@@ -411,7 +415,8 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
     else if (dmrginp.specificpdm().size()==2)
       npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[1]);
     else abort();
-    p3out << "\t\t\t NPDM sweep time " << timerX.elapsedwalltime() << " " << timerX.elapsedcputime() << endl;
+    ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime();
+    p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl;
     return;
   }
 
@@ -427,14 +432,16 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
 				dmrginp.npdm_generate() = true;
         SweepGenblock::do_one(sweepParams, false, !direction, false, 0, state, stateB); //this will generate the cd operators                               
 				dmrginp.npdm_generate() = false;
-        p3out << "\t\t\t NPDM SweepGenblock time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+        ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+        p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl;
         dmrginp.set_fullrestart() = false;
 
         sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy;
         Timer timerX;
         npdm_driver->clear();
         npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,stateB);
-        p3out << "\t\t\t NPDM sweep time " << timerX.elapsedwalltime() << " " << timerX.elapsedcputime() << endl;
+        ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime();
+        p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl;
        }
        }
     }
@@ -446,7 +453,8 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
     sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy;
     SweepGenblock::do_one(sweepParams, false, !direction, false, 0, -1, -1); //this will generate the cd operators                               
 		dmrginp.npdm_generate() = false;
-    p3out << "\t\t\t NPDM SweepGenblock time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+    ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+    p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl;
     dmrginp.set_fullrestart() = false;
 
 
@@ -460,7 +468,8 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
           Timer timerX;
           npdm_driver->clear();
           npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,stateB);
-          p3out << "\t\t\t NPDM sweep time " << timerX.elapsedwalltime() << " " << timerX.elapsedcputime() << endl;
+          ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime();
+          p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl;
           }
       }
 
@@ -472,7 +481,8 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
           Timer timerX;
           npdm_driver->clear();
           npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,state);
-          p3out << "\t\t\t NPDM sweep time " << timerX.elapsedwalltime() << " " << timerX.elapsedcputime() << endl;
+          ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime();
+          p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl;
         } 
         else{
           if (npdm_order == NPDM_ONEPDM) SweepOnepdm::do_one(sweepParams, false, direction, false, 0, state);     
@@ -516,13 +526,15 @@ void npdm(NpdmOrder npdm_order, bool restartpdm, bool transitionpdm)
         SweepGenblock::do_one(sweepParams, false, !direction, false, 0, state, stateB); //this will generate the cd operators
 				dmrginp.npdm_generate() = false;
         dmrginp.set_fullrestart() = false;
-        p3out << "\t\t\t NPDM SweepGenblock time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+        ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+        p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl;
 
         Timer timerX;
         npdm_driver->clear();
         sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy;
         npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0, state,stateB);
-        p3out << "\t\t\t NPDM sweep time " << timerX.elapsedwalltime() << " " << timerX.elapsedcputime() << endl;
+        ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime();
+        p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl;
       }
     }
   }

modules/npdm/npdm_driver.C

@@ -719,7 +719,8 @@ void Npdm_driver::compute_npdm_elements(std::vector<Wavefunction> & wavefunction
     loop_over_operator_patterns_store( npdm_patterns, npdm_expectations, big );
 #ifndef SERIAL
   world.barrier();
-  p3out << " prepare intermediate time"<< timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+  ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+  p3out << " prepare intermediate time"<< ewall << " " << ecpu << endl;
 #endif
   }
 
@@ -772,15 +773,9 @@ void Npdm_driver::compute_npdm_elements(std::vector<Wavefunction> & wavefunction
   } else {
     reduce(world, write_intermediate_time, std::plus<double>(), 0);
   }
-  if (mpigetrank() == 0) {
-    double sum, sum2;
-    reduce(world, timer.elapsedwalltime(), sum2, std::plus<double>(), 0);
-    reduce(world, timer.elapsedcputime(), sum, std::plus<double>(), 0);
-    p3out << "NPDM compute elements time " << sum << " "<< sum2 << endl;
-  } else {
-    reduce(world, timer.elapsedwalltime(), std::plus<double>(), 0);
-    reduce(world, timer.elapsedcputime(), std::plus<double>(), 0);
-  }
+
+  ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+  p3out << "NPDM compute elements time " << ewall << " "<< ecpu << endl;
 #else
   p3out << "NPDM compute elements time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
 #endif

modules/npdm/onepdm_container.C

@@ -52,7 +52,8 @@ void Onepdm_container::save_npdms(const int& i, const int& j)
 #ifndef SERIAL
   world.barrier();
 #endif
-  p3out << "1PDM save full array time " << timer.elapsedwalltime() << " " << timer.elapsedcputime() << endl;
+  ecpu = timer.elapsedcputime();ewall=timer.elapsedwalltime();
+  p3out << "1PDM save full array time " << ewall << " " << ecpu << endl;
 }
 
 //-----------------------------------------------------------------------------------------------------------------------------------------------------------