-
Notifications
You must be signed in to change notification settings - Fork 369
/
jellium_hf_state.py
99 lines (80 loc) · 4.3 KB
/
jellium_hf_state.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""This module constructs the uniform electron gas' Hartree-Fock state."""
# TODO: Move this to linalg after circular import issue has been fixed.
import numpy
from openfermion.ops.operators import FermionOperator
from openfermion.transforms.opconversions import normal_ordered
from openfermion.transforms.repconversions import inverse_fourier_transform
from openfermion.utils import count_qubits
from openfermion.hamiltonians.jellium import plane_wave_kinetic
def lowest_single_particle_energy_states(hamiltonian, n_states):
"""Find the lowest single-particle states of the given Hamiltonian.
Args:
hamiltonian (FermionOperator)
n_states (int): Number of lowest energy states to give."""
# Enumerate the single-particle states.
n_single_particle_states = count_qubits(hamiltonian)
# Compute the energies for each of the single-particle states.
single_particle_energies = numpy.zeros(n_single_particle_states, dtype=float)
for i in range(n_single_particle_states):
single_particle_energies[i] = hamiltonian.terms.get(((i, 1), (i, 0)), 0.0)
# Find the n_states lowest states.
occupied_states = single_particle_energies.argsort()[:n_states]
return list(occupied_states)
def hartree_fock_state_jellium(grid, n_electrons, spinless=True, plane_wave=False):
"""Give the Hartree-Fock state of jellium.
Args:
grid (Grid): The discretization to use.
n_electrons (int): Number of electrons in the system.
spinless (bool): Whether to use the spinless model or not.
plane_wave (bool): Whether to return the Hartree-Fock state in
the plane wave (True) or dual basis (False).
Notes:
The jellium model is built up by filling the lowest-energy
single-particle states in the plane-wave Hamiltonian until
n_electrons states are filled.
"""
# Get the jellium Hamiltonian in the plane wave basis.
# For determining the Hartree-Fock state in the PW basis, only the
# kinetic energy terms matter.
hamiltonian = plane_wave_kinetic(grid, spinless=spinless)
hamiltonian = normal_ordered(hamiltonian)
hamiltonian.compress()
# The number of occupied single-particle states is the number of electrons.
# Those states with the lowest single-particle energies are occupied first.
occupied_states = lowest_single_particle_energy_states(hamiltonian, n_electrons)
occupied_states = numpy.array(occupied_states)
if plane_wave:
# In the plane wave basis the HF state is a single determinant.
hartree_fock_state_index = numpy.sum(2**occupied_states)
hartree_fock_state = numpy.zeros(2 ** count_qubits(hamiltonian), dtype=complex)
hartree_fock_state[hartree_fock_state_index] = 1.0
else:
# Inverse Fourier transform the creation operators for the state to get
# to the dual basis state, then use that to get the dual basis state.
hartree_fock_state_creation_operator = FermionOperator.identity()
for state in occupied_states[::-1]:
hartree_fock_state_creation_operator *= FermionOperator(((int(state), 1),))
dual_basis_hf_creation_operator = inverse_fourier_transform(
hartree_fock_state_creation_operator, grid, spinless
)
dual_basis_hf_creation = normal_ordered(dual_basis_hf_creation_operator)
# Initialize the HF state.
hartree_fock_state = numpy.zeros(2 ** count_qubits(hamiltonian), dtype=complex)
# Populate the elements of the HF state in the dual basis.
for term in dual_basis_hf_creation.terms:
index = 0
for operator in term:
index += 2 ** operator[0]
hartree_fock_state[index] = dual_basis_hf_creation.terms[term]
return hartree_fock_state