1 body reduced density matrix estimator produces incorrect results

[jptowns]

Collecting the "dm1b" estimator on a VMC calculation of FeO with a single determinant PBE orbital basis with no Jastrow factor should produce a diagonal density matrix. This system contains 22 electrons per spin channel, so a visualization of the 1rdm for, say, the "up" channel should show a strongly diagonal structure for the first 22 states. Yet, for QMCPACK built on 04 Sept. (#1046) it does not:

The salient input block is in the specification of the hamiltonian:

<?xml version="1.0"?>                                                                                                  
<qmcsystem>                                                                                                            
  <hamiltonian name="h0" type="generic" target="e">                                                                    
    <pairpot name="IonIon"   type="coulomb" source="i" target="i"/>                                                    
    <pairpot name="ElecElec" type="coulomb" source="e" target="e"/>                                                    
    <pairpot type="pseudo" name="PseudoPot" source="i" wavefunction="psi0" format="xml">                               
      <pseudo elementType="Fe" href="../wfns-and-pseudos/Fe.opt.xml"/>                                                 
      <pseudo elementType="O"  href="../wfns-and-pseudos/O.opt.xml"/>                                                  
    </pairpot>                                                                                                         
    <estimator type="dm1b" name="1rdms">                                                                               
      <parameter name="basis"            >  dm_basis     </parameter>                                                  
      <parameter name="points"           >  10           </parameter>                                                  
      <parameter name="samples"          >  10           </parameter>                                                  
      <parameter name="warmup"           >  1            </parameter>                                                  
      <parameter name="timestep"         >  0.25         </parameter>                                                  
    </estimator>                                                                                                       
  </hamiltonian>                                                                                                       
</qmcsystem>                                                                                                           

And the VMC calculation was specified by:

  <qmc method="vmc" move="pbyp" checkpoint="0">                                                                        
    <parameter name="walkers"       >       8  </parameter>                                                            
    <parameter name="timestep"      >  0.2000  </parameter>                                                            
    <parameter name="useDrift"      >      no  </parameter>                                                            
    <parameter name="steps"         >      64  </parameter>                                                            
    <parameter name="blocks"        >      64  </parameter>                                                            
    <parameter name="nonlocalmoves" >     yes  </parameter>                                                            
  </qmc>                                                                                                               

What is strange is that for commit #1021 the same input blocks produce something more expected:

In fact, these effects persist regardless of the presence of a Jastrow, and also persist at the DMC level.
Minimal working example included.

feo-mwe-1rdm.gz

[prckent]

Is the trace correct?

[jptowns]

Whoa! The trace is way off!
The trace of the rdms for the first figure (As of commit #1046) are about 11 and 0, for up and down spin channels, respectively. Should be ~22 and ~22.

The trace for the second figure (as of commit #1021) are 21.9/21.9 respectively. This is what is expected.

[jptowns]

Oh, also I know that there are currently some issues regarding OMP support so all calculations were performed on a desktop workstation with multiple MPI ranks each with 1 thread.

[jtkrogel]

Interesting. There aren't many PR's between #1021 and #1046 that could affect this. Have you narrowed down to the culprit PR (via bisection search, etc)?

[jptowns]

Have not done that. Will see what I can do.

[prckent]

Any updates on this?

[prckent]

Is this fixed by #1509 or are there still problems to be resolved?

[jtkrogel]

I will look into this soon. I suspect it may be a real vs. complex issue. The real orbitals are not orthogonal in general and might lead to the pattern in the 1rdm seen above.

[prckent]

Could we measure using a distinct complex-valued orbital set or an orthogonal real-valued set even with the real build?

[jtkrogel]

If the set of orbitals is orthogonal within QMCPACK, then everything should work fine. I would prefer to require orthogonal orbitals within QMCPACK overall.

[Paul-St-Young]

Is there any update on this issue?

I still get 11 and 0 as the traces for the 1RDMs of up and down electrons using Josh's example above.

                    QMCPACK 3.9.9

       (c) Copyright 2003-  QMCPACK developers

                    Please cite:
 J. Kim et al. J. Phys. Cond. Mat. 30 195901 (2018)
      https://doi.org/10.1088/1361-648X/aab9c3

  Git branch: develop
  Last git commit: 02efd65977f97fc8ca54340640701553741445de
  Last git commit date: Tue Sep 22 08:04:13 2020 -0500

[prckent]

I recall: (1) it is still an issue (2) there are no spin polarized 1RDM tests (statistical or deterministic) (3) no one complained enough to get this fixed sooner. I recall that we do check the trace these days. Clearly a plumbing issue.

[jtkrogel]

@Paul-St-Young have you tried with the complex build of QMCPACK?

[jtkrogel]

Also, would you report the output with <parameter name="check_overlap" > yes </parameter>?

[Paul-St-Young]

@jtkrogel this run was performed using complex build.
The overlap matrix is within 1e-4 of being the identity.
feo-dm_basis-ovlp.zip

[jtkrogel]

But VMC w/o a Jastrow does not yield the identity in the same channel as the selected basis?

[Paul-St-Young]

correct.

[jtkrogel]

Would you mind uploading an ascii file containing the mean of the sampled density matrices?

[Paul-St-Young]

Here are the means of the up and down DMs
feo-mew-1rdm-dats.zip
and the outputs from plot_rdm.py for a quick peak (first up then down).
The down DM is certainly wrong. Feels like a memory address scramble to me.

[jtkrogel]

Interesting. Possibly not a coincidence that the diagonal dominance of the sub-blocks fits a 2x2 pattern (i.e. falling in exactly the last half of the 30 orbitals).

I agree this looks like an indexing problem.

[jtkrogel]

Do you have time to bisect between #1021 and #1046?

[Paul-St-Young]

@jtkrogel I tested out the commit from #1021. The RDMs still look the same as above so I did not attach them.

Last git commit: 0b5eab33ec8fa981d44fc51fcb7652fc289d6a4c
Last git commit date: Tue Aug 21 09:34:14 2018 -0400
Last git commit subject: Merge pull request #1021 from markdewing/cusp_structure

I did find that the problem goes away if the evaluator is switched from "loop" to "matrix".
I guess maybe default it to "matrix" and remove the "loop" evaluator?
1021_matrix.zip

[Paul-St-Young]

In DensityMatrix1B.cpp::evaluate_loop, int ij = nindex + s * basis_size2; is inside the ns loop. This looks suspicious.

[Paul-St-Young]

NVM, moving int ij = nindex + s * basis_size2; outside creates seg. faults...
Is the "loop" evaluator worth fixing?

[jtkrogel]

If loop is still the default, it should be switched to matrix. Loop is the original implementation, and matrix is much more efficient (it's what I've always used).

In light of future ports (i.e. gpu/offload compatibility) perhaps not taking the full matrix approach, I think loop should remain, but an abort added noting the current bug. We can make the decision later if fixing the loop implementation makes sense.

Thanks for determining the difference.

[aannabe]

Sharing some results to isolate the problem. Using recent builds #4060, 8-atom Si supercell with 16-up + 16-down electrons, the following is the trace summary for various dm1b parameters:

evaluator   integrator      WaveFunc      Trace_Up
======================================================
matrix      uniform         Slater        8.001(2)
matrix      uniform_grid    Slater        7.99998(1)
matrix      density         Slater        15.995(9)
loop        uniform         Slater        7.999(2)
loop        uniform_grid    Slater        7.99999(3)
loop        density         Slater        15.999(9)

matrix      density         Slater-Jast   15.71(1)

Correct Value:                            16

Basis: 32 PBE  orbitals

Both loop and matrix evaluator traces are correct for only the density integrator. Others produce only half of the electron number. Below are the visualizations for matrix evaluator using various integrators:

matrix with uniform:

matrix with uniform_grid:

matrix with density:

Including Jastrow introduces the "smearing" of the diagonal as expected.
matrix with density and Jastrows:

Therefore, the density integrator seems to be working correctly at the moment. The uniform and uniform_grid 1RDMs look different even though their traces are the same.

Runs:
si_bulk.zip

[jtkrogel]

What is printed when you run with check_overlap on?

[aannabe]

Here is the plot of the overlap matrix. It's mostly diagonal, but there are off-diagonal terms as large as 0.6. These are PBE orbitals.