QMCPACK / QMCPACK/qmcpack

Add Ghost-Atom Support in QMCPACK/NEXUS.

Open
#5,570 4 comments 0 reactions 0 assignees View on GitHub
enhancement stale
Dominant language
C++
Stars
403
Forks
154
Avg merge
1d 12h
Merged PRs (30d)
82

Description

**Feature request:** Add ghost‐atom support to PyscfToQmcpack.py, convert4qmc, and NEXUS so that QMCPACK can perform counterpoise BSSE corrections.

**Rationale:** Basis‐set superposition error (BSSE) can dominate binding‐energy errors. Counterpoise corrections with modest basis sets (e.g. cc-pVTZ, aug-cc-pVDZ) approach the complete‐basis‐set limit (Nakano *et al.*, JCTC 2025) [DOI:10.1021/acs.jctc.4c01631]. Primer about BSSE [here](http://vergil.chemistry.gatech.edu/notes/cp.pdf).

**Current workaround:** I apply an ad-hoc patch (ccECP only on real atoms) to reproduce Nakano’s results. Details below; would love native ghost-atom handling instead.

## 1. To make PySCF’s ghost atom system work with QMCPACK.
The only script that needs modification is `PyscfToQmcpack.py`, so that it can correctly interpret ghost atoms from the PySCF output.

Since ghost atoms are defined simply by a prefix (e.g. `"X-"`), they can be handled by adding entries to the `IonName` dictionary with an atomic number of zero.

In `PyscfToQmcpack.py`, add:
```
IonName.update({f"X-{el}": 0 for el in ["H","C"]})
```

## 2. Generate & fix QMCPACK XML
Assuming no errors appear from the `savetoqmcpack` function in `PyscfToQmcpack.py`, the next step is to use `convert4qmc` to generate QMCPACK files.

### 2.1. Edit `.structure.xml`
The `.structure.xml` file needs some editing. `convert4qmc` generates an empty `` element in the `particleset` group:

```

0
0
0

```

When the conversion routine encounters a species with atomic number zero, it creates a blank `` entry. We’ll replace that placeholder with explicit ghost-atom groups.

For example:
```

0
0
0

0
0
0

```

The `attrib name="position"` attribute prints the correct atomic positions, including ghost atoms, but `attrib name="ionid"` leaves out all ghost atoms, so we need to add them again.

For example:
```


X-C X-C X-H X-H X-H X-H H C H H H

```
Verify that the number of `` entries matches `` rows to avoid misalignment.

### 2.2. Update `` in `.qmc.in`
Modify the ECP part to ensure the `` blocks in the Hamiltonian section match ccECP specs and that ghost species carry no ccECP entry.

## 3. Add ghost-atom Jastrow terms
I found that inserting these correlation blocks for each ghost species makes the Jastrow optimization work as expected.

For example:
```


```

This setup yields a variance-to-energy ratio for the ghost-atom calculation that is essentially the same as for the system without ghost atoms.

## 4. Run DMC (no changes)
No further edits are required for the DMC step. You can use the exact same .qmc.in settings and workflow you’d apply to the standard system—just point to your ghost-augmented wavefunction and structure files. In our tests, the only difference is that the Hamiltonian and Jastrow blocks already include the ghost species, but all other DMC parameters (time step, number of walkers, equilibration steps, etc.) remain identical.

## 5. Extend NEXUS for ghost system
As you can see from sections 1 to 3, the calculation process is similar to what NEXUS produces for a standard system. Therefore, it was possible for NEXUS to automate everything by editing the NEXUS Python library. Here, I will explain which library I edited and why.

### 5.1 Recognize ghosts in `periodic_table.py`
The first error is NEXUS does not recognize ghost atom like `X-H` or `X-C`.

on class `SimpleElement`:
```diff
class SimpleElement(DevBase):
def __init__(self):
...
+ self.ghost = False
...
+ def copy(self):
+ new = deepcopy(self)
+ return new
```

on class `PeriodicTable`:
```diff
class PeriodicTable(DevBase):
+ GHOST_PREFIX = 'X-' # or '*' or any convention you like

element_set=set([
'Ac','Al','Am','Sb','Ar','As','At','Ba','Bk','Be','Bi','B' ,'Br',
'Cd','Ca','Cf','C' ,'Ce','Cs','Cl','Cr','Co','Cu','Cm','Dy','Es',
...
isotopes[symbol] = elem_isotopes
#end for
self.isotopes = isotopes

+ # … build self.simple_elements and self.elements exactly as before …
+ # but remove the hacky `self[elem.symbol] = element` and instead:
+ self._base = self.elements # real symbols → Element objects

#end def __init__
...
#end class PeriodicTable
```
With this change, there’s no need to edit the huge dictionary every time you add a new ghost atom; the ghost atom prefix can also be changed.

### 5.2 Propagate ghosts in `physical_system.py`
This time `PhysicalSystem` shows this error:
```
PhysicalSystem error:
particle X-H is unknown
exiting.
```

In `physical_system.py`, edit the `PhysicalSystem` class:
```diff
def add_particles(self, **particle_counts):
pc = self.particle_collection # all known particles
plist = []
for name,count in particle_counts.items():
+ # detect ghost-prefix (e.g. "X-" or your chosen prefix)
+ ghost = False
+ if name.startswith("X-"):
+ ghost = True
+ lookup_name = name[2:] # strip "X-"
+ else:
+ lookup_name = name
particle = pc.get_particle(lookup_name)
if particle is None:
self.error('particle {0} is unknown'.format(name))
else:
particle = particle.copy()
+ setattr(particle, "ghost", ghost) # mark it ghost for downstream code
#end if
particle.set_count(count)
plist.append(particle)
#end for
self.particles.add_particles(plist)
self.update()
#end def add_particles
```

This will:
1. Checks if your incoming name starts with `"X-"`.
2. Strips off the prefix before calling `get_particle`.
4. Copies the real particle, then sets `particle.ghost=True`.

### 5.3 Handle ghost species in `qmcpack_input.py`
Now the error is more specific:
```
You said:
Traceback (most recent call last):
File "/home/s2420025/Volumes/storage/hdd4/reza/06ghost/12test/02dmc_monoA/script.py", line 129, in
optJ2 = generate_qmcpack(
File "/home/s2420025/app/qmcpack-3.17.1/nexus/lib/qmcpack.py", line 1591, in generate_qmcpack
sim_args.input = generate_qmcpack_input(**inp_args)
File "/home/s2420025/app/qmcpack-3.17.1/nexus/lib/qmcpack_input.py", line 7375, in generate_qmcpack_input
inp = generate_basic_input(**kwargs)
File "/home/s2420025/app/qmcpack-3.17.1/nexus/lib/qmcpack_input.py", line 7586, in generate_basic_input
particlesets = generate_particlesets(
File "/home/s2420025/app/qmcpack-3.17.1/nexus/lib/qmcpack_input.py", line 4566, in generate_particlesets
ion = ions[ion_spec]
File "/home/s2420025/app/qmcpack-3.17.1/nexus/lib/generic.py", line 158, in __getitem__
return self.__dict__[name]
KeyError: 'X-C'
```

I proceeded to edit generate_particlesets in qmcpack_input.py:
```diff
def generate_particlesets(electrons = 'e',
ions = 'ion0',
up = 'u',
...
groups = []
+ ghost = False
for ion_spec in ion_species:
- ion = ions[ion_spec]
+ # strip off ghost-prefix if present
+ ghost = False
+ key = ion_spec
+ if key.startswith("X-"):
+ ghost = True
+ key = key[2:]
+ # now look up the real ion
+ base_ion = ions[key]
+ ion = base_ion.copy() # clone so we don’t mutate the original
+ ion.ghost = ghost # mark it ghost if needed
gpos = pos[elem==ion.name]
g = group(
name = ion.name,
...
```
After this, when `ion_spec` is `"X-C"`, it will strip `"X-"`, look up `"C"`, copy it, flag it as a ghost, and carry on.

## 6. Remaining manual step.
With these changes in place, NEXUS will recognize and propagate ghost-atom species through almost every stage of the QMC workflow, from species creation and input file generation to Jastrow setup and the full DMC run. However, Even after these changes, NEXUS stop right after `convert4qmc` and you must still manually edit .structure.xml and .qmc.in. Automating that final patch would close the loop.

## 7. Example calculation: BSSE evaluation on CH4⋯⋯C2H4 with increasing basis set

**Figure 1** shows the convergence of the DMC binding energy for the CH4⋯⋯C2H4 van der Waals complex toward the CBS limit (computed with aug-cc-pV6Z) as we increase the augmented correlation-consistent basis set. For each basis set, aug-cc-pVDZ and aug-cc-pVTZ, we plot both the raw DMC bonding energy (no CP) and the counterpoise corrected result (CP). Error bars are one-sigma QMC uncertainties. The plain binding energy at aug-cc-pVDZ underbinds by approximately 0.03 eV; applying the CP correction recovers nearly the CBS value (within approximately 0.01 eV). By aug-cc-pVTZ the uncorrected error is already around 0.02 eV.

Raw calculation data for this example is [here](https://www.dropbox.com/scl/fi/w34w8zqnbtjqwfedwwnqg/QMCPACK_BSSE.tar.gz?rlkey=0hmen0qnck1u9i6lr9q64sxnw&st=gucbchcp&dl=0).

Image

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.