Evaluating desired estimators without evaluating kinetic energy/local energy

[Ch-Y Wang]

Dear QMCpack developers,

Since I need to evaluate observables such as spin densities with huge amount of samples in order to reduce statistical error (particularly for smaller bin size), I need to surpass the default local energy estimator to reduce the computational cost. The local energy esitmator in this case is not important any more because I already have the optimized Jastrow parameters from previous VMC runs. 

I currently run into two difficulties: the mainestimator is by default the local energy estimator, and it seems the code still doing local energy calculations even I switched it to SpinDensityNew. For the hamiltonian part, I tried to set "default = "no" , but QMCpack  reduced to run with this setting.

Please find my inputs in the attachment. Thank you very much.

Best wishes,

Chung-Yu 

[Paul R. C. Kent]

Dear Chung Yu -

I have been trying to find a complete solution to this that will not
require any changes to the current code. The challenge is to find a
way to run VMC with only estimators being evaluated, and no
evaluations of the energy/Hamiltonian. Depending on the cost of the
estimators this will significantly speed up the run. As far as I know,
no one else has tried to run this way in a long time, but happily it
does turn out to be possible.

After some experimentation I have a working solution, but it does use
a workaround: while it looks like the both the modern batched VMC
driver and the legacy VMC driver should be able to run and not
evaluate any Hamiltonian, this is not the case due to a
bug. Additionally, the kinetic energy is implicitly added and this
needs to be removed. The workaround is to leave a trivial constant
part of the Hamiltonian present.

My experiments were based on running the Nexus workflow example
qmcpack/nexus/examples/qmcpack/rsqmc_quantum_espresso/01_diamond_dft_vmc/diamond_lda_vmc_batched_estimators.py
and, changing it to only compute the spindensity, and after running
successfully, trying to cut down further the input *.in.xml for the
VMC run. Similar experiments with the legacy drivers obtained the same
result. Note that the estimators are pulled out of the Hamiltonian XML
block in the the batched drivers for exactly this kind of calculation.

QMCPACK enforces that there is a <hamiltonian></hamiltonian> present
in the input and will abort if not. You can remove the non-local
potentials and electron-ion and electron-electron Coulomb terms
leaving an empty <hamiltonian> section. In my testing, this works
correctly. However, the kinetic energy (only) is still evaluated. As
described in
https://qmcpack.readthedocs.io/en/develop/hamiltonianobservable.html#the-hamiltonian
, this is added implicitly. Adding the attribute default="no" to the
Hamiltonian block avoids this implicit addition of the kinetic energy
term, but results in a crash:

```
  Adding psi0 TrialWaveFunction to the pool

 Hamiltonian and observables
 ---------------------------
  Name: h0

std::bad_alloc
Fatal Error. Aborting at Unhandled Exception
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.
```

After some experimentation it became clear that we have a bug when the
Hamiltonian is completely empty.

e.g. If we we have a typical
Hamiltonian block and put default="no",
then the Kinetic energy is correctly not evaluated and there are no
problems with the other terms or the estimators

```

      <hamiltonian name="h0" type="generic" target="e" default="no">
         <pairpot type="coulomb" name="ElecElec" source="e" target="e"/>
         <pairpot type="coulomb" name="IonIon" source="ion0" target="ion0"/>
         <pairpot type="pseudo" name="PseudoPot" source="ion0" wavefunction="psi0" format="xml">
            <pseudo elementType="C" href="C.BFD.xml"/>
         </pairpot>
         <pairpot type="MPC" name="MPC" source="e" target="e" ecut="60.0" physical="no"/>
      </hamiltonian>
      <estimators>
         <estimator type="spindensity" name="SpinDensity">
            <parameter name="grid">
               8 8 8
            </parameter>
         </estimator>
      </estimators>
   </qmcsystem>
   <qmc method="vmc" move="pbyp">
      <parameter name="walkers_per_rank"    >    4               </parameter>
      <parameter name="warmupSteps"         >    50              </parameter>
      <parameter name="blocks"              >    5               </parameter>
      <parameter name="steps"               >    5               </parameter>
      <parameter name="subSteps"            >    3               </parameter>
      <parameter name="timestep"            >    0.3             </parameter>
      <parameter name="useDrift"            >    no              </parameter>
   </qmc>
```
Then the run completes with the kinetic energy correctly not
         evaluated:

```
cat vmc.s000.scalar.dat
#   index    LocalEnergy         LocalEnergy_sq      LocalPotential      LocalECP            NonLocalECP         ElecElec            IonIon              MPC                 BlockWeight         BlockCPU            AcceptRatio        
         0   -8.2142615254e+01    6.8187626704e+03   -5.7847877227e+01   -2.4294738027e+01    4.2461791159e+00   -1.0991377781e+01   -5.1102678561e+01   -1.0744261965e+01    2.0000000000e+01    2.0098958000e-02    5.3854166667e-01
         1   -8.0901448557e+01    6.5735790756e+03   -5.6605002424e+01   -2.4296446133e+01    5.2327450054e+00   -1.0735068868e+01   -5.1102678561e+01   -1.0432414531e+01    2.0000000000e+01    2.2338083000e-02    5.6718750000e-01
         2   -7.9958736709e+01    6.4628972478e+03   -5.7604630935e+01   -2.2354105775e+01    3.8725697883e+00   -1.0374522161e+01   -5.1102678561e+01   -9.9840382727e+00    2.0000000000e+01    2.3065833000e-02    5.5260416667e-01
         3   -7.8547504500e+01    6.2084813982e+03   -5.8538660003e+01   -2.0008844498e+01    4.2877830942e+00   -1.1723764535e+01   -5.1102678561e+01   -1.1439376486e+01    2.0000000000e+01    2.2834750000e-02    5.5260416667e-01
         4   -8.0925701383e+01    6.7662941289e+03   -5.3757205627e+01   -2.7168495756e+01    5.5693895063e+00   -8.2239165714e+00   -5.1102678561e+01   -7.8843445237e+00    2.0000000000e+01    2.3424250000e-02    5.4947916667e-01
```
The workaround is to have just the ion-ion energy evaluated:
```

      <hamiltonian name="h0" type="generic" target="e" default="no">
         <pairpot type="coulomb" name="IonIon" source="ion0" target="ion0"/>
      </hamiltonian>
      <estimators>
         <estimator type="spindensity" name="SpinDensity">
            <parameter name="grid">
               8 8 8
            </parameter>
         </estimator>
      </estimators>
   </qmcsystem>
   <qmc method="vmc" move="pbyp">
      <parameter name="walkers_per_rank"    >    4               </parameter>
      <parameter name="warmupSteps"         >    50              </parameter>
      <parameter name="blocks"              >    5               </parameter>
      <parameter name="steps"               >    5               </parameter>
      <parameter name="subSteps"            >    3               </parameter>
      <parameter name="timestep"            >    0.3             </parameter>
      <parameter name="useDrift"            >    no              </parameter>
   </qmc>
</simulation>
```
which results in a scalar file like:
```
$ cat vmc.s000.scalar.dat
#   index    LocalEnergy         LocalEnergy_sq      LocalPotential      IonIon              BlockWeight         BlockCPU            AcceptRatio        
         0   -5.1102678561e+01    2.6114837562e+03    0.0000000000e+00   -5.1102678561e+01    2.0000000000e+01    1.5569625000e-02    5.6354166667e-01
         1   -5.1102678561e+01    2.6114837562e+03    0.0000000000e+00   -5.1102678561e+01    2.0000000000e+01    1.6128417000e-02    5.5312500000e-01
         2   -5.1102678561e+01    2.6114837562e+03    0.0000000000e+00   -5.1102678561e+01    2.0000000000e+01    1.5904459000e-02    5.7291666667e-01
         3   -5.1102678561e+01    2.6114837562e+03    0.0000000000e+00   -5.1102678561e+01    2.0000000000e+01    1.5990958000e-02    5.5989583333e-01
         4   -5.1102678561e+01    2.6114837562e+03    0.0000000000e+00   -5.1102678561e+01    2.0000000000e+01    1.5196500000e-02    5.5052083333e-01
```
while the density is in the stat.h5 file:
```
$ h5ls -r vmc.s000.stat.h5
/                        Group
/IonIon                  Group
/IonIon/value            Dataset {5/Inf, 1}
/LocalEnergy             Group
/LocalEnergy/value       Dataset {5/Inf, 1}
/LocalEnergy_sq          Group
/LocalEnergy_sq/value    Dataset {5/Inf, 1}
/LocalPotential          Group
/LocalPotential/value    Dataset {5/Inf, 1}
/SpinDensity             Group
/SpinDensity/d           Group
/SpinDensity/d/value     Dataset {5/Inf, 512}
/SpinDensity/u           Group
/SpinDensity/u/value     Dataset {5/Inf, 512}
```

The "empty Hamiltonian" bug probably occurs due to a simple underlying assumption that there
is at least one term in the Hamiltonian somewhere.

Hopefully this solves your problem.

-- Paul