frozen-core in dual plane wave basis/'pseudo-potential'/valence-electrons only

[hrwernerer]

I have a question regarding the dual plane wave basis.
In my understanding I can create a Fermionic hamiltonian describing the electronic structure in a periodic material, on a discretized grid, with a combination of the Hamiltonian constructors dual_basis_jellium_model() and adding dual_basis_external_potential() for the electron-nucleic interactions. (The latter function has a look-up table to simply count the positive charge number for the given molecular geometry nuclei.)

As I am mostly interested in semiconductors/metals, I would next like to apply the valence electron or 'frozen core' method, as a parallel to the workflow in molecular chemistry where we may as an approximation freeze the lowest-energy orbitals and for example only consider valence electrons. For that, we would normally use freeze_orbitals() on the Fermionic Hamiltonianin in f.ex. a slater-type orbital basis set.

In the case of dual plane wave basis however, it does not make as much sense (I believe), because the Fermionic modes do not refer to spatial 'orbitals' but rather plane waves, right?

I think what I'm looking for is a way to incorporate ab-initio "pseudo-potentials" for the nuclei, and only consider the valence electrons in the Hamiltonian operator, with the quantization into a dual plane wave basis.
Is such functionality currently featured in OpenFermion, is there a plan to add it, or could you point me in the right direction to contribute such a thing?

[babbush]

Indeed, if you call freeze_orbitals() on the plane wave or dual basis Hamiltonians it will just freeze the lowest energy plane wave or dual basis functions. But what you would really want to do is something like freeze the lowest energy molecular orbitals. For that you would basically need to be working in a molecular orbital basis. However you are also correct that a viable alternative to this which would apply in the case of the plane wave or dual basis sets is to pseudopotentials. Currently, we do not have code for using pseudopotentials in OpenFermion but it is a reasonable thing to want to have. Basically, the external potential operator would need to be computed differently if working with pseudopotentials. Exactly how it would be computed depends on the pseudopotential. Note that for many classes of pseudopotentials, including e.g., non-local pseudopotentials, the external potential operator will no longer be diagonal in the dual basis, but that I don't expect that would be much of a problem. So the first step would be to derive an expression for the integrals defining the external potential. I.e.,
U = \sum_{pq} U_{pq} a^\dagger_p a_q
where U_{pq} = \int dr \phi_p(r) U(r) \phi_q(r)^*
where \phi_p(r) is the spatial function of plane wave with mode q and U(r) is the pseudopotential of interest. Specifically, if we are dealing with atom centered pseudopotentials then U(r) is a sum of the pseudopotentials centered on each atom, e.g., U(r) = \sum_{j} W_{Z_j} (r - R_j) where W_Z(r) is the atom centered pseudopotential centered at the origin for the atom with charge Z, Z_j is the charge of atom j and R_j is the position of atom j.