Issue with 1-RDM calculation using VMC for molecules

[yuanliu1]

Describe the bug I am trying to compute the 1-RDM of a H2 molecule under STO-3G basis from VMC runs, by setting the Hamiltonian block in the input file as the following. The vmc run ends successfully and produced the .stat.h5 file with a "DensityMatrices" of size (200,2,2) with all entries being 0. Two questions: 1) The 1-RDM under STO-3G basis should be of size 2x2, and I don't understand why the stat.h5 file contains a (200,2,2) object; 2) Why all entries of the stat.h5 file are zero?

Is something wrong with my input? I am attaching all the input files below.

  <hamiltonian name="h0" type="generic" target="e">
  <estimator type="dm1b" name="DensityMatrices">
  <parameter name="basis"        >  updet downdet  </parameter>
  <parameter name="evaluator"    >  matrix        </parameter>
  <parameter name="integrator"   >  uniform_grid  </parameter>
  <parameter name="points"       >  10             </parameter>
  <parameter name="scale"        >  1.0           </parameter>
  <parameter name="center"       >  0 0 0         </parameter>
  </estimator>
    <pairpot name="ElecElec" type="coulomb" source="e" target="e" physical="true"/>
    <pairpot name="IonIon" type="coulomb" source="ion0" target="ion0"/>
    <pairpot name="IonElec" type="coulomb" source="ion0" target="e"/>
  </hamiltonian>

  <qmc method="vmc" move="pbyp" checkpoint="-1">
    <estimator name="LocalEnergy" hdf5="no"/>
    <parameter name="warmupSteps">100</parameter>
    <parameter name="blocks">200</parameter>
    <parameter name="steps">50</parameter>
    <parameter name="substeps">8</parameter>
    <parameter name="timestep">0.5</parameter>
    <parameter name="usedrift">no</parameter>
    <!--Sample count should match targetwalker count for DMC. Will be obtained from all nodes.-->
    <parameter name="samples">16000</parameter>
  </qmc>
</simulation>

To Reproduce Steps to reproduce the behavior: I am running qmcpack as a docker container (williamfgc/qmcpack-ci:ubuntu20-openmpi) as in https://qmcpack.readthedocs.io/en/develop/running_docker.html. It passes all unit tests.

Run qmcpack input.xml

Expected behavior Will generate stat.h5 file which is supposed to contain the 1-RDM.

Additional context Only using J2 factor in this case for simplicity. The PySCF script used to generate the wave function file is also attached. h2_test.zip

[jtkrogel]

In answer to your questions:

1) The size is (200,2,2) because there are 200 statistical samples of the matrix (equal to number of blocks). You will need to take the mean over the blocks (after excluding an equlibration period) prior to e.g. diagonalizing the matrix to get natural orbitals and occupancies. Resampling methods like jackknife or bootstrap can be used to obtain errorbars on the occupancies.

2. The integration is performed on a uniform grid, per your input. In this case, you will need to define a spatial cell to perform the additional integration over. An old example for an oxygen atom is below. Another (less tested) option you might try is to set integrator to "density" and optionally provide values for samples/warmup/timestep in the density matrix input (the defaults should be OK though).

Density matrix section in manual: https://qmcpack.readthedocs.io/en/develop/hamiltonianobservable.html#one-body-density-matrix

Old example input for an oxygen atom using open boundary conditions, but with a cell defined for the integration grid:

<simulation>
  <project id="qmc" series="0">
    <application name="qmcapp" role="molecu" class="serial" version="0.2">
      QMC of F2 molecule, Roos Triple Zeta basis
    </application>
  </project>

  <qmcsystem>
    <simulationcell name="global">
      <parameter name="lattice" units="bohr">
        16.0   0.0   0.0
         0.0  16.0   0.0
         0.0   0.0  16.0
      </parameter>
      <parameter name="bconds">
        n n n
      </parameter>
    </simulationcell>
  </qmcsystem>

  <include href="O.Gaussian-G2.ptcl.xml"/>
  <include href="hf.Gaussian-G2.cusp.sposet.xml"/>

  <hamiltonian name="h0" type="generic" default="multi" target="e">
    <pairpot name="ElecElec" type="coulomb" source="e" target="e"/>
    <pairpot name="Coulomb" type="coulomb" source="ion0" target="ion0"/>
    <pairpot name="Coulomb" type="coulomb" source="ion0" target="e"/>

    <estimator type="dm1b" name="DensityMatrices">
      <parameter name="energy_matrix"       >    yes                   </parameter>
      <parameter name="integrator"          >    uniform_grid          </parameter>
      <parameter name="points"              >    10                    </parameter>
      <parameter name="scale"               >    1.0                   </parameter>
      <parameter name="basis"               >    dm_basis              </parameter>
      <parameter name="normalized"          >    yes                   </parameter>
      <parameter name="evaluator"           >    matrix                </parameter>
      <parameter name="center"              >    0 0 0                 </parameter>
      <parameter name="check_overlap"       >    no                    </parameter>
      <parameter name="rstats"              >    no                    </parameter>
    </estimator>

  </hamiltonian>
  <qmc method="vmc" move="pbyp" target="e">
    <parameter name="walkers">       1 </parameter>
    <parameter name="warmupsteps"> 200 </parameter>
    <parameter name="blocks">      800 </parameter>
    <parameter name="steps">        10 </parameter>
    <parameter name="timestep">    0.5 </parameter>
    <parameter name="samples">     960 </parameter>
  </qmc>
  <qmc method="dmc" move="pbyp" target="e">
    <parameter name="warmupsteps"> 20 </parameter>
    <parameter name="blocks">      10 </parameter>
    <parameter name="steps">        5 </parameter>
    <parameter name="timestep">   0.04 </parameter>
  </qmc>
  <qmc method="dmc" move="pbyp" target="e">
    <parameter name="warmupsteps"> 20  </parameter>
    <parameter name="blocks">     400  </parameter>
    <parameter name="steps">        5  </parameter>
    <parameter name="timestep">   0.01 </parameter>
  </qmc>
</simulation>

[yuanliu1]

Hi Jaron,

Thanks for your answer. Setting a simulation cell solved the problem - now I get some non-zero elements in the stat.h5 file (the first few lines) but those numbers does not agree with the true 1-RDM (below) computed from PySCF for my Hamiltonian: from stat.h5 output GROUP "/" { GROUP "DensityMatrices" { GROUP "number_matrix" { GROUP "d" { DATASET "value" { DATATYPE H5T_IEEE_F64LE DATASPACE SIMPLE { ( 200, 2, 2 ) / ( H5S_UNLIMITED, 2, 2 ) } DATA { (0,0,0): 0.000402533, 0.000402533, (0,1,0): 0.000402533, 0.000402533, (1,0,0): 0.000400865, 0.000400865, (1,1,0): 0.000400865, 0.000400865, (2,0,0): 0.000401147, 0.000401147, (2,1,0): 0.000401147, 0.000401147, (3,0,0): 0.000404192, 0.000404192, (3,1,0): 0.000404192, 0.000404192, From PySCF: 1-RDM = [[0.70873275 0.58665908] [0.58665908 0.70873275]].

I think maybe because I passed the wrong basis in the following `

updet downdet ` whereas the oxygen example above uses some sposet generated elsewhere ` dm_basis `. So I was trying to build a sposet using the following code, but the code errors out with `Fatal Error. Aborting at BasisSetBase::createSPOSet(cur,input_info) has not been implemented` seems to be suggesting some functionaliy of creating the sposet was not implemented. Is there a way to solve this? ` 16.0 0.0 0.0 0.0 16.0 0.0 0.0 0.0 16.0 n n n 0.7473652629 0.3732680729 0.1330985646 -0.02509993107 -0.1330034768 -0.2216508189 -0.2916191937 -0.3525409865 -0.192285904 -0.03446824124 `

[jtkrogel]

Can you post a zipped archive with files to reproduce your problem?