Cusp correction only works if atoms of same type are consecutive

[zenandrea]

Describe the bug
Cusp correction (in the CASINO way, for all-electron calculations, with a GTO basis set) only works if atoms of same type are consecutive. For instance, if I give the input of a molecule in the following way:
H -0.923961000000000 1.231956000000000 -1.684781230000000
C 0.000000000000000 0.667176000000000 -1.684781230000000
C 0.000000000000000 -0.667176000000000 -1.684781230000000
H 0.923961000000000 1.231956000000000 -1.684781230000000
H 0.923961000000000 -1.231956000000000 -1.684781230000000
H -0.923961000000000 -1.231956000000000 -1.684781230000000
The cusp correction gives totally wrong energies, while if I input the molecule in the following way:
C 0.000000000000000 0.667176000000000 -1.684781230000000
C 0.000000000000000 -0.667176000000000 -1.684781230000000
H -0.923961000000000 1.231956000000000 -1.684781230000000
H 0.923961000000000 1.231956000000000 -1.684781230000000
H 0.923961000000000 -1.231956000000000 -1.684781230000000
H -0.923961000000000 -1.231956000000000 -1.684781230000000
the cusp correction works properly.

To Reproduce
Steps to reproduce the behavior:

1. release version or git commit hash being built: observed in the last 3.17.1 version, as well 3.14 and 3.15.
2. cmake command
3. full program/test invocation command
4. additional steps

Expected behavior
A clear and concise description of what you expected to happen.

I expect the cusp correction should work independently on the order of atoms, or the code should stop with an error message when the atom order is not consistent with the cusp correction.

System:

• system name: both on my laptop (mac os x) and on Archer2 cluster
• modules loaded [e.g. output of module list]
• other systems where this is reproducible [e.g. "my laptop", "none"]

Additional context
Add any other context about the problem here.

[prckent]

Thanks for the report. I guess it has always been this way. Should be simple to make independent of atom order or (worst case) abort.

[markdewing]

See #4145 (and linked issues and fixes) for previous encounters with cusp correction and atom reordering.

[prckent]

Thanks Mark. Had forgotten our previous encounters here.

[ye-luo]

I guess you have input files for both cases. Could you provide them? So we can reproduce the exact issue.

[zenandrea]

Dear @prckent and @markdewing, it took me some time and effort to discover the source of the unreasonable energies I was initially obtaining. If you did know about the issue with cusp-correction since July 2022, why not adding a check that stops the calculation when there is cusp-correction and the atoms are not ordered?
As a user, it really saves a lot of time if the run just stops when the input has something that prevents it from working properly.

Anyway, while preparing a reproducible example for @ye-luo, I actually discovered that the problem does not always appear with “disordered atoms”, sometimes it does, and sometimes it doesn’t.
In particular, I was observing the problem when I was using nexus, while I do not see the problem if I use convert4qmc (and qmcpack version 3.17.1, because with version 3.14 I still observe the problem).

This is an example of a system test:
run a pyscf calculation as the one copied below (named scf.py), followed by

python ./scf.py | tee scf.out

convert4qmc -addCusp -nojastrow  -orbitals scf.h5
cp scf.wfnoj.xml _scf.wfnoj.xml
cat _scf.wfnoj.xml | sed 's/..\/updet/updet/g' | sed 's/..\/downdet/downdet/g' > scf.wfnoj.xml

qmcpack scf.qmc.in-wfnoj.xml

grep converged scf.out
qmca -q ev *scalar.dat

This is the python script (scf.py) to generate the orbitals on scf.h5:

#! /usr/bin/env python

import numpy as np
from numpy import array
from pyscf import df, scf, dft
### generated system text ###
from pyscf import gto as gto_loc
mol = gto_loc.Mole()
mol.verbose  = 5
mol.atom     = '''
  C   -1.061709204   1.297140572   0.292060003
  O   -0.358161116   2.270458613   0.531812668
  O   -0.589303516   0.094917758   0.003788813
  H    0.404435659   0.127722621   0.018411838
  C   -2.558427798   1.342549823   0.296257320
  H   -2.895997978   2.347464002   0.518316340
  H   -2.932889278   1.022390451  -0.672995551
  H   -2.937211960   0.644910433   1.039557084
'''
mol.basis    = 'ccpvdz'
mol.unit     = 'A'
mol.charge   = 0
mol.spin     = 0
mol.symmetry = False
mol.build()
mf = scf.RHF(mol)
e_scf = mf.kernel()
### generated conversion text ###
from PyscfToQmcpack import savetoqmcpack
savetoqmcpack(mol,mf,'scf')

and an alternative calculation can be run with consecutive atoms

mol.atom     = '''
  O   -0.358161116   2.270458613   0.531812668
  O   -0.589303516   0.094917758   0.003788813
  C   -1.061709204   1.297140572   0.292060003
  C   -2.558427798   1.342549823   0.296257320
  H    0.404435659   0.127722621   0.018411838
  H   -2.895997978   2.347464002   0.518316340
  H   -2.932889278   1.022390451  -0.672995551
  H   -2.937211960   0.644910433   1.039557084
'''

As the scf.py is an HF calculation, we can compare the energy from pyscf and qmcpack (for a vmc without Jastrow), expecting only a little difference because of the cusp-correction.

For the input below I get:

$ grep converged scf.out
converged SCF energy = -227.829505206272
$ qmca -q ev *scalar.dat
                            LocalEnergy               Variance           ratio
scf  series 0  -227.653796 +/- 0.216931   51.139856 +/- 3.746079   0.2246

If I use the version 3.14 of qmcpack I get with the “diserdered” atoms:

scf  series 0  -90.268001 +/- 15.531538   1101263.449672 +/- 302640.165916   12199.9317

while obtaining the same energy from pyscf.
However, in v. 3.17.1 (but also v. 3.15 on my laptop) I get the same numbers from both calculations.

With nexus I instead obtain a totally wrong energy also with the latest version.
So, I looked at the differences in the qmcpack input when the cuspInfo files are generated and I noticed that the issue might be related with the definition of the particleset, which is the following (this is another molecule I tested)

$ diff TEST_cuspcorr*/cusp.in.xml
23,46c23,45
<       <particleset name="ion0">
<          <group name="H" size="4" mass="1837.3622193391143">
<             <parameter name="charge"              >    1                     </parameter>
<             <parameter name="valence"             >    1                     </parameter>
<             <parameter name="atomicnumber"        >    1                     </parameter>
<             <parameter name="mass"                >    1837.3622193391143            </parameter>
<             <attrib name="position" datatype="posArray" condition="0">
<                     -1.74603325        2.32805945       -3.18377512
<                      1.74603325        2.32805945       -3.18377512
<                      1.74603325       -2.32805945       -3.18377512
<                     -1.74603325       -2.32805945       -3.18377512
<             </attrib>
<          </group>
<          <group name="C" size="2" mass="21894.71359057295">
<             <parameter name="charge"              >    6                     </parameter>
<             <parameter name="valence"             >    6                     </parameter>
<             <parameter name="atomicnumber"        >    6                     </parameter>
<             <parameter name="mass"                >    21894.71359057295            </parameter>
<             <attrib name="position" datatype="posArray" condition="0">
<                      0.00000000        1.26077992       -3.18377512
<                      0.00000000       -1.26077992       -3.18377512
<             </attrib>
<          </group>
<       </particleset>
---
>   <particleset name="ion0" size="6">
>     <group name="H">
>       <parameter name="charge">1</parameter>
>       <parameter name="valence">1</parameter>
>       <parameter name="atomicnumber">1</parameter>
>     </group>
>     <group name="C">
>       <parameter name="charge">6</parameter>
>       <parameter name="valence">6</parameter>
>       <parameter name="atomicnumber">6</parameter>
>     </group>
>     <attrib name="position" datatype="posArray">
>  -1.7460332398e+00  2.3280594375e+00 -3.1837751045e+00
>   0.0000000000e+00  1.2607799169e+00 -3.1837751045e+00
>   0.0000000000e+00 -1.2607799169e+00 -3.1837751045e+00
>   1.7460332398e+00  2.3280594375e+00 -3.1837751045e+00
>   1.7460332398e+00 -2.3280594375e+00 -3.1837751045e+00
>  -1.7460332398e+00 -2.3280594375e+00 -3.1837751045e+00
> </attrib>
>     <attrib name="ionid" datatype="stringArray">
>  H C C H H H
> </attrib>
>   </particleset>

The particleset below works correctly, the one above doesn’t, and I think that comes from the fact that in the input above the overall order of atoms on the underlying orbital file (scf.h5, or ../pySCF/c4q.orbs.h5) is not given.

By the way: convert4qmc generated an input such that the cusp-correction files are on the parent directory (../) instead of the present working directory... which I think it is not a good idea.

[prckent]

| If you did know about the issue with cusp-correction since July 2022...
According to the linked issues, the problem was thought to be solved back then.

[zenandrea]

@prckent I'm sorry, I thought it was still to be fixed. So, this one is a "new" bug.

[ye-luo]

In the $ diff TEST_cuspcorr*/cusp.in.xml, both the above and below styles are supported.
The above one first puts 4 H and then 2 C and it mismatches your orbs.h5.
The below one H C C H H H specifies the exact placement. That is what you expected.

In brief pyscf->convert4qmc->qmcpack works fine. pyscf->nexus->qmcpack fails?
I think either the info of atoms specified in the nexus input wasn't right or nexus doesn't extract the atoms from orbs.h5 correctly. Could you provide your nexus script?

By the way: convert4qmc generated an input such that the cusp-correction files are on the parent directory (../) instead of the present working directory... which I think it is not a good idea.

I agree.

[zenandrea]

Ok, the nexus script below yields correct results for H2O.xyz equal to:

3

O  0.000000  0.000000  0.000000
H  0.000000  0.757160  0.586260
H  0.000000  0.757160 -0.586260`

and totally wrong for:

3

H  0.000000  0.757160  0.586260
O  0.000000  0.000000  0.000000
H  0.000000  0.757160 -0.586260`

#! /usr/bin/env python3

from nexus import settings,job,run_project,obj
from nexus import generate_physical_system
from nexus import generate_pyscf
from nexus import generate_convert4qmc
from nexus import generate_cusp_correction
from nexus import generate_qmcpack

settings(
    results    = '',
    sleep      = 3,
    machine    = 'ws16',
    )

system = generate_physical_system(
    structure = 'H2O.xyz',
    )

# perform Hartree-Fock
scf = generate_pyscf(
    identifier = 'scf',               # log output goes to scf.out
    path       = 'H2O/hf',            # directory to run in
    job        = job(serial=True),    # pyscf must run serially
    template   = './scf_template.py', # pyscf template file
    system     = system,
    mole       = obj(                 # used to make Mole() inputs
        verbose  = 5,
        basis    = 'ccpvtz',
        symmetry = True,
        ),
    save_qmc   = True,                # save wfn data for qmcpack
    )

# convert orbitals to QMCPACK format
c4q = generate_convert4qmc(
    identifier   = 'c4q',
    path         = 'H2O/hf',
    job          = job(cores=1),
    no_jastrow   = True,
    hdf5         = True,              # use hdf5 format
    dependencies = (scf,'orbitals'),
    )

# calculate cusp correction
cc = generate_cusp_correction(
    identifier   = 'cusp',
    path         = 'H2O/cuspcorr',
    job          = job(cores=16,threads=16),
    system       = system,
    dependencies = (c4q,'orbitals'),
    )

# collect dependencies relating to orbitals
orbdeps = [(c4q,'particles'), # pyscf changes particle positions
           (c4q,'orbitals' ),
           (cc ,'cuspcorr' )]

# vmc without Jastrow
vmcJ0 = generate_qmcpack(
    #block             = True,
    identifier        = 'vmc',
    path              = 'H2O/vmcJ0',
    job               = job(cores=8),
    system            = system,
    jastrows = [],
    qmc               = 'vmc',        # quartic variance optimization
    dependencies      = orbdeps,
    )

run_project()