/
commutator_diagonal_coulomb_operator.py
312 lines (258 loc) · 13.3 KB
/
commutator_diagonal_coulomb_operator.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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
# 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.
"""Faster commutators for two-body operators with diagonal Coulomb terms."""
import warnings
from openfermion.ops.operators import FermionOperator
from openfermion.transforms.opconversions import term_reordering
def commutator_ordered_diagonal_coulomb_with_two_body_operator(
operator_a, operator_b, prior_terms=None
):
"""Compute the commutator of two-body operators provided that both are
normal-ordered and that the first only has diagonal Coulomb interactions.
Args:
operator_a: The first FermionOperator argument of the commutator.
All terms must be normal-ordered, and furthermore either hopping
operators (i^ j) or diagonal Coulomb operators (i^ i or i^ j^ i j).
operator_b: The second FermionOperator argument of the commutator.
operator_b can be any arbitrary two-body operator.
prior_terms (optional): The initial FermionOperator to add to.
Returns:
The commutator, or the commutator added to prior_terms if provided.
Notes:
The function could be readily extended to the case of arbitrary
two-body operator_a given that operator_b has the desired form;
however, the extra check slows it down without desirable added utility.
"""
if prior_terms is None:
prior_terms = FermionOperator.zero()
for term_a in operator_a.terms:
coeff_a = operator_a.terms[term_a]
for term_b in operator_b.terms:
coeff_b = operator_b.terms[term_b]
coefficient = coeff_a * coeff_b
# If term_a == term_b the terms commute, nothing to add.
if term_a == term_b or not term_a or not term_b:
continue
# Case 1: both operators are two-body, operator_a is i^ j^ i j.
if (
len(term_a) == len(term_b) == 4
and term_a[0][0] == term_a[2][0]
and term_a[1][0] == term_a[3][0]
):
_commutator_two_body_diagonal_with_two_body(
term_a, term_b, coefficient, prior_terms
)
# Case 2: commutator of a 1-body and a 2-body operator
elif (len(term_b) == 4 and len(term_a) == 2) or (len(term_a) == 4 and len(term_b) == 2):
_commutator_one_body_with_two_body(term_a, term_b, coefficient, prior_terms)
# Case 3: both terms are one-body operators (both length 2)
elif len(term_a) == 2 and len(term_b) == 2:
_commutator_one_body_with_one_body(term_a, term_b, coefficient, prior_terms)
# Final case (case 4): violation of the input promise. Still
# compute the commutator, but warn the user.
else:
warnings.warn(
'Defaulted to standard commutator evaluation ' 'due to an out-of-spec operator.'
)
additional = FermionOperator.zero()
additional.terms[term_a + term_b] = coefficient
additional.terms[term_b + term_a] = -coefficient
additional = term_reordering.normal_ordered(additional)
prior_terms += additional
return prior_terms
def _commutator_one_body_with_one_body(
one_body_action_a, one_body_action_b, coefficient, prior_terms
):
"""Compute the commutator of two one-body operators specified by actions.
Args:
one_body_action_a, one_body_action_b (tuple): single terms of one-body
FermionOperators (i^ j or i^ i).
coefficient (complex float): coefficient of the commutator.
prior_terms (FermionOperator): prior terms to add the commutator to.
"""
# In the case that both the creation and annihilation operators of the
# two actions pair, two new terms must be added.
if (
one_body_action_a[0][0] == one_body_action_b[1][0]
and one_body_action_b[0][0] == one_body_action_a[1][0]
):
new_one_body_action_a = ((one_body_action_a[0][0], 1), (one_body_action_a[0][0], 0))
new_one_body_action_b = ((one_body_action_b[0][0], 1), (one_body_action_b[0][0], 0))
prior_terms.terms[new_one_body_action_a] = (
prior_terms.terms.get(new_one_body_action_a, 0.0) + coefficient
)
prior_terms.terms[new_one_body_action_b] = (
prior_terms.terms.get(new_one_body_action_b, 0.0) - coefficient
)
# A single pairing causes the mixed action a[0]^ b[1] to be added
elif one_body_action_a[1][0] == one_body_action_b[0][0]:
action_ab = ((one_body_action_a[0][0], 1), (one_body_action_b[1][0], 0))
prior_terms.terms[action_ab] = prior_terms.terms.get(action_ab, 0.0) + coefficient
# The other single pairing adds the mixed action b[0]^ a[1]
elif one_body_action_a[0][0] == one_body_action_b[1][0]:
action_ba = ((one_body_action_b[0][0], 1), (one_body_action_a[1][0], 0))
prior_terms.terms[action_ba] = prior_terms.terms.get(action_ba, 0.0) - coefficient
def _commutator_one_body_with_two_body(action_a, action_b, coefficient, prior_terms):
"""Compute commutator of action-specified one- and two-body operators.
Args:
action_a, action_b (tuple): single terms, one one-body and the other
two-body, from normal-ordered FermionOperators. It does not matter
which is one- or two-body so long as only one of each appears.
coefficient (complex float): coefficient of the commutator.
prior_terms (FermionOperator): prior terms to add the commutator to.
"""
# Determine which action is 1-body and which is 2-body.
# Label the creation and annihilation parts of the two terms.
if len(action_a) == 4 and len(action_b) == 2:
one_body_create = action_b[0][0]
one_body_annihilate = action_b[1][0]
two_body_create = (action_a[0][0], action_a[1][0])
two_body_annihilate = (action_a[2][0], action_a[3][0])
new_action = list(action_a)
# Flip coefficient because we reversed the commutator's arguments.
coefficient *= -1
else:
one_body_create = action_a[0][0]
one_body_annihilate = action_a[1][0]
two_body_create = (action_b[0][0], action_b[1][0])
two_body_annihilate = (action_b[2][0], action_b[3][0])
new_action = list(action_b)
# If both terms are composed of number operators, they commute.
if one_body_create == one_body_annihilate and (two_body_create == two_body_annihilate):
return
# If the one-body annihilation is in the two-body creation parts
if one_body_annihilate in two_body_create:
new_coeff = coefficient
new_inner_action = list(new_action)
# Determine which creation part(s) of the one-body action to use
if one_body_annihilate == two_body_create[0]:
new_inner_action[0] = (one_body_create, 1)
elif one_body_annihilate == two_body_create[1]:
new_inner_action[1] = (one_body_create, 1)
# Normal order if necessary
if new_inner_action[0][0] < new_inner_action[1][0]:
new_inner_action[0], new_inner_action[1] = (new_inner_action[1], new_inner_action[0])
new_coeff *= -1
# Add the resulting term.
if new_inner_action[0][0] > new_inner_action[1][0]:
prior_terms.terms[tuple(new_inner_action)] = (
prior_terms.terms.get(tuple(new_inner_action), 0.0) + new_coeff
)
# If the one-body creation is in the two-body annihilation parts
if one_body_create in two_body_annihilate:
new_coeff = -coefficient
# Determine which annihilation part(s) of the one-body action to sub in
if one_body_create == two_body_annihilate[0]:
new_action[2] = (one_body_annihilate, 0)
elif one_body_create == two_body_annihilate[1]:
new_action[3] = (one_body_annihilate, 0)
# Normal order if necessary
if new_action[2][0] < new_action[3][0]:
new_action[2], new_action[3] = new_action[3], new_action[2]
new_coeff *= -1
# Add the resulting term.
if new_action[2][0] > new_action[3][0]:
prior_terms.terms[tuple(new_action)] = (
prior_terms.terms.get(tuple(new_action), 0.0) + new_coeff
)
def _commutator_two_body_diagonal_with_two_body(
diagonal_coulomb_action, arbitrary_two_body_action, coefficient, prior_terms
):
"""Compute the commutator of two two-body operators specified by actions.
Args:
diagonal_coulomb_action (tuple): single term of a diagonal Coulomb
FermionOperator (i^ j^ i j). Must be in normal-ordered form,
i.e. i > j.
arbitrary_two_body_action (tuple): arbitrary single term of a two-body
FermionOperator, in normal-ordered form, i.e. i^ j^ k l with
i > j, k > l.
coefficient (complex float): coefficient of the commutator.
prior_terms (FermionOperator): prior terms to add the commutator to.
Notes:
The function could be readily extended to the case of reversed input
arguments (where diagonal_coulomb_action is the arbitrary one, and
arbitrary_two_body_action is from a diagonal Coulomb FermionOperator);
however, the extra check slows it down without significantly increased
utility.
"""
# Identify creation and annihilation parts of arbitrary_two_body_action.
arb_2bdy_create = (arbitrary_two_body_action[0][0], arbitrary_two_body_action[1][0])
arb_2bdy_annihilate = (arbitrary_two_body_action[2][0], arbitrary_two_body_action[3][0])
# The first two sub-cases cover when the creations and annihilations of
# diagonal_coulomb_action and arbitrary_two_body_action totally pair up.
if (
diagonal_coulomb_action[2][0] == arbitrary_two_body_action[0][0]
and diagonal_coulomb_action[3][0] == arbitrary_two_body_action[1][0]
):
prior_terms.terms[arbitrary_two_body_action] = (
prior_terms.terms.get(arbitrary_two_body_action, 0.0) - coefficient
)
elif (
diagonal_coulomb_action[0][0] == arbitrary_two_body_action[2][0]
and diagonal_coulomb_action[1][0] == arbitrary_two_body_action[3][0]
):
prior_terms.terms[arbitrary_two_body_action] = (
prior_terms.terms.get(arbitrary_two_body_action, 0.0) + coefficient
)
# Exactly one of diagonal_coulomb_action's creations matches one of
# arbitrary_two_body_action's annihilations.
elif diagonal_coulomb_action[0][0] in arb_2bdy_annihilate:
# Nothing gets added if there's an unbalanced double creation.
if diagonal_coulomb_action[1][0] in arb_2bdy_create or (
diagonal_coulomb_action[0][0] in arb_2bdy_create
):
return
_add_three_body_term(
arbitrary_two_body_action, coefficient, diagonal_coulomb_action[1][0], prior_terms
)
elif diagonal_coulomb_action[1][0] in arb_2bdy_annihilate:
# Nothing gets added if there's an unbalanced double creation.
if diagonal_coulomb_action[0][0] in arb_2bdy_create or (
diagonal_coulomb_action[1][0] in arb_2bdy_create
):
return
_add_three_body_term(
arbitrary_two_body_action, coefficient, diagonal_coulomb_action[0][0], prior_terms
)
elif diagonal_coulomb_action[0][0] in arb_2bdy_create:
_add_three_body_term(
arbitrary_two_body_action, -coefficient, diagonal_coulomb_action[1][0], prior_terms
)
elif diagonal_coulomb_action[1][0] in arb_2bdy_create:
_add_three_body_term(
arbitrary_two_body_action, -coefficient, diagonal_coulomb_action[0][0], prior_terms
)
def _add_three_body_term(two_body_action, coefficient, mode, prior_terms):
new_action = list(two_body_action)
# Insert creation and annihilation operators into the two-body action.
new_action.insert(0, (mode, 1))
new_action.insert(3, (mode, 0))
# Normal order the creation operators of the new action.
# Each exchange in the action flips the sign of the coefficient.
if new_action[0][0] < new_action[1][0]:
new_action[0], new_action[1] = new_action[1], new_action[0]
coefficient *= -1
if new_action[1][0] < new_action[2][0]:
new_action[1], new_action[2] = new_action[2], new_action[1]
coefficient *= -1
# Normal order the annihilation operators of the new action.
# Each exchange in the action flips the sign of the coefficient.
if new_action[3][0] < new_action[4][0]:
new_action[3], new_action[4] = new_action[4], new_action[3]
coefficient *= -1
if new_action[4][0] < new_action[5][0]:
new_action[4], new_action[5] = new_action[5], new_action[4]
coefficient *= -1
# Add the new normal-ordered term to the prior terms.
prior_terms.terms[tuple(new_action)] = (
prior_terms.terms.get(tuple(new_action), 0.0) + coefficient
)