generator.py 11.8 KB
Newer Older
Davison's avatar
Davison committed
1
2
import tensorflow as tf
# tf.enable_v2_behavior()
3
4
from tensornetwork import *
import numpy as np
Davison's avatar
Davison committed
5
6
from oop_imsrg.hamiltonian import *
from oop_imsrg.occupation_tensors import *
7

8
9
10
11
12
13
14
15
class Generator(object):
    """Parent class for organization purposes. Ideally, all Generator
    classes should inherit from this class. In this way, AssertionErrors
    can be handled in a general way."""
    def calc_eta():
        print("Function that calculates the generator")

class WegnerGenerator(Generator):
Davison's avatar
Davison committed
16
17
    """Calculate Wegner's generator for a normal ordered Hamiltonian.
       Truncate at two-body interactions."""
18
19

    def __init__(self, h, occ_t):
20
21
22
23
24
25
        """Class constructor. Instantiate WegnerGenerator object.

        Arguments:

        h -- Hamiltonian object (must be normal-ordered)
        occ_t -- OccupationTensor object"""
Davison's avatar
Davison committed
26

27
        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
        assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"

        self.f = h.f
        self.G = h.G

        self._holes = h.holes
        self._particles = h.particles

        self._occA = occ_t.occA
        self._occB = occ_t.occB
        self._occC = occ_t.occC
        self._occD = occ_t.occD

    @property
    def f(self):
43
44
45
        """Returns:

        f -- one-body tensor elements (initialized by Hamiltonian object)"""
46
47
48
49
        return self._f

    @property
    def G(self):
50
51
        """Returns:

Davison's avatar
Davison committed
52
        f -- two-body tensor elements (initialized by Hamiltonian object)"""
53
54
55
56
        return self._G

    @f.setter
    def f(self, f):
57
        """Sets the one-body tensor."""
58
59
60
61
        self._f = f

    @G.setter
    def G(self, G):
62
        """Sets the two-body tensor."""
63
        self._G = G
Davison's avatar
Davison committed
64

65
    # @classmethod
Davison's avatar
Davison committed
66
    def decouple_OD(self):
Davison's avatar
Davison committed
67
68
        """Decouple the off-/diagonal elements from each other in
        the one- and two-body tensors. This procedure is outlined in
69
70
71
72
73
74
75
76
77
        An Advanced Course in Computation Nuclear Physics, Ch.10.

        Returns:

        (fd, -- diagonal part of f
         fod, -- off-diagonal part of f
         Gd, -- diagonal part of G
         God) -- off-diagonal part of G"""

78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
        f = self.f
        G = self.G
        holes = self._holes
        particles = self._particles

        # - Decouple off-diagonal 1B and 2B pieces
        fod = np.zeros(f.shape)
        fod[np.ix_(particles, holes)] += f[np.ix_(particles, holes)]
        fod[np.ix_(holes, particles)] += f[np.ix_(holes, particles)]
        fd = f - fod

        God = np.zeros(G.shape)
        God[np.ix_(particles, particles, holes, holes)] += G[np.ix_(particles, particles, holes, holes)]
        God[np.ix_(holes, holes, particles, particles)] += G[np.ix_(holes, holes, particles, particles)]
        Gd = G - God

        return (fd, fod, Gd, God)

96
    # @classmethod
97
    def calc_eta(self):
Davison's avatar
Davison committed
98
        """Calculate the generator. The terms are defined in An
99
100
101
102
103
104
105
106
        Advanced Course in Computation Nuclear Physics, Ch.10.
        See also dx.doi.org/10.1016/j.physrep.2015.12.007

        Returns:

        (eta1B, -- one-body generator
         eta2B) -- two-body generator"""

107
108
109
110
111
        partition = self.decouple_OD()
        fd = partition[0]
        fod = partition[1]
        Gd = partition[2]
        God = partition[3]
112
113
114
115
116

        holes = self._holes
        particles = self._particles

        occA = self._occA
117
        occB = self._occB
118
119
120
121
122
123
124
        occC = self._occC

        # - Calculate 1B generator
        # first term
        sum1_1b_1 = ncon([fd, fod], [(-1, 0), (0, -2)]).numpy()
        sum1_1b_2 = np.transpose(sum1_1b_1)
        sum1_1b = sum1_1b_1 - sum1_1b_2
Davison's avatar
Davison committed
125

126
127
128
129
130
        # second term
        sum2_1b_1 = ncon([fd, God], [(0, 1), (1, -1, 0, -2)]).numpy()
        sum2_1b_2 = ncon([fod, Gd], [(0, 1), (1, -1, 0, -2)]).numpy()
        sum2_1b_3 = sum2_1b_1 - sum2_1b_2
        sum2_1b = ncon([occA, sum2_1b_3],[(-1, -2, 0, 1), (0,1)]).numpy()
Davison's avatar
Davison committed
131

132
133
134
135
136
        # third term
        sum3_1b_1 = ncon([occC, God], [(-1, -2, -3, 0, 1, 2), (0, 1, 2, -4)]).numpy()
        sum3_1b_2 = ncon([Gd, sum3_1b_1], [(2, -1, 0, 1), (0, 1, 2, -2)]).numpy()
        sum3_1b_3 = np.transpose(sum3_1b_2)
        sum3_1b = sum3_1b_2 - sum3_1b_3
Davison's avatar
Davison committed
137

138
139
140
141
142
143
144
145
146
        eta1B = sum1_1b + sum2_1b + 0.5*sum3_1b

        # - Calculate 2B generator
        # first term (P_ij piece)
        sum1_2b_1 = ncon([fd, God], [(-1, 0), (0, -2, -3, -4)]).numpy()
        sum1_2b_2 = ncon([fod, Gd], [(-1, 0), (0, -2, -3, -4)]).numpy()
        sum1_2b_3 = sum1_2b_1 - sum1_2b_2
        sum1_2b_4 = np.transpose(sum1_2b_3, [1, 0, 2, 3])
        sum1_2b_5 = sum1_2b_3 - sum1_2b_4
Davison's avatar
Davison committed
147

148
149
150
151
152
153
        # first term (P_kl piece)
        sum1_2b_6 = ncon([fd, God], [(0, -3), (-1, -2, 0, -4)]).numpy()
        sum1_2b_7 = ncon([fod, Gd], [(0, -3), (-1, -2, 0, -4)]).numpy()
        sum1_2b_8 = sum1_2b_6 - sum1_2b_7
        sum1_2b_9 = np.transpose(sum1_2b_8, [0, 1, 3, 2])
        sum1_2b_10 = sum1_2b_8 - sum1_2b_9
Davison's avatar
Davison committed
154

155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
        sum1_2b = sum1_2b_5 - sum1_2b_10

        # second term
        sum2_2b_1 = ncon([occB, God], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
        sum2_2b_2 = ncon([occB,  Gd], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
        sum2_2b_3 = ncon([Gd,  sum2_2b_1], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
        sum2_2b_4 = ncon([God, sum2_2b_2], [(-1, -2, 0, 1), (0, 1, -3, -4)]).numpy()
        sum2_2b = sum2_2b_3 - sum2_2b_4

        # third term
        sum3_2b_1 = ncon([Gd, God], [(0, -1, 1, -3), (1, -2, 0, -4)]).numpy()
        sum3_2b_2 = np.transpose(sum3_2b_1, [1, 0, 2, 3])
        sum3_2b_3 = np.transpose(sum3_2b_1, [0, 1, 3, 2])
        sum3_2b_4 = np.transpose(sum3_2b_1, [1, 0, 3, 2])
        sum3_2b_5 = sum3_2b_1 - sum3_2b_2 - sum3_2b_3 + sum3_2b_4
        sum3_2b = ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)]).numpy()

        eta2B = sum1_2b + 0.5*sum2_2b + sum3_2b

174
175
        return (eta1B, eta2B)

Davison's avatar
Davison committed
176
class WegnerGenerator3B(WegnerGenerator):
177

Davison's avatar
Davison committed
178
179
    """Calculate Wegner's generator for a normal ordered Hamiltonian.
       Truncate at three-body interactions. Inherits from WegnerGenerator."""
180

Davison's avatar
Davison committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
    def __init__(self, h, occ_t):
        """Class constructor. Instantiate WegnerGenerator object.

        Arguments:

        h -- Hamiltonian object (must be normal-ordered)
        occ_t -- OccupationTensor object"""

        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
        assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"

        self.f = h.f
        self.G = h.G
        self.W = np.zeros(h.n_sp_states*np.ones(6,dtype=np.int64))

        self._holes = h.holes
        self._particles = h.particles

        self._occA = occ_t.occA
Jacob Davison's avatar
Jacob Davison committed
200
        self._occA2 = occ_t.occA2
Davison's avatar
Davison committed
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
        self._occB = occ_t.occB
        self._occC = occ_t.occC
        self._occD = occ_t.occD
        self._occE = occ_t.occE
        self._occF = occ_t.occF
        self._occG = occ_t.occG
        self._occH = occ_t.occH
        self._occI = occ_t.occI
        self._occJ = occ_t.occJ

    @property
    def W(self):
        """Returns:

        W -- three-body tensor elements (initialized by Hamiltonian object)"""
        return self._W

    @W.setter
    def W(self, W):
        """Sets the three-body tensor."""
        self._W = W

223
    # @classmethod
Davison's avatar
Davison committed
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
    def decouple_OD(self):
        """Inherits from WegnerGenerator.

        Decouple the off-/diagonal elements from each other in
        the one- and two-body tensors. This procedure is outlined in
        An Advanced Course in Computation Nuclear Physics, Ch.10.

        Returns:

        (fd, -- diagonal part of f
         fod, -- off-diagonal part of f
         Gd, -- diagonal part of G
         God, -- off-diagonal part of G
         Wd, -- diagonal part of W
         Wod) -- off-diagonal part of W"""
239

240
241
242
243
244
        partition = super().decouple_OD()
        fd = partition[0]
        fod = partition[1]
        Gd = partition[2]
        God = partition[3]
Davison's avatar
Davison committed
245

Davison's avatar
Davison committed
246
        W = self.W
247

Davison's avatar
Davison committed
248
249
250
251
252
253
254
255
256
257
        holes = self._holes
        particles = self._particles

        Wod = np.zeros(W.shape)
        Wod[np.ix_(particles, particles, particles, holes, holes, holes)] += W[np.ix_(particles, particles, particles, holes, holes, holes)]
        Wod[np.ix_(holes, holes, holes, particles, particles, particles)] += W[np.ix_(holes, holes, holes, particles, particles, particles)]
        Wd = W - Wod

        return(fd, fod, Gd, God, Wd, Wod)

258
    # @classmethod
Davison's avatar
Davison committed
259
260
261
    def calc_eta(self):
        """Inherits from WegnerGenerator.

262
263
        Calculate the generator. See dx.doi.org/10.1016/j.physrep.2015.12.007,
        Appendix B, for three-body flow equations.
Davison's avatar
Davison committed
264
265
266
267
268
269
270

        Returns:

        (eta1B, -- one-body generator
         eta2B, -- two-body generator
         eta3B) -- three-body generator"""

271
272
273
274
275
276
277
        partition = self.decouple_OD()
        fd = partition[0]
        fod = partition[1]
        Gd = partition[2]
        God = partition[3]
        Wd = partition[4]
        Wod = partition[5]
Davison's avatar
Davison committed
278

279
        eta1B, eta2B = super().calc_eta()
Davison's avatar
Davison committed
280
281
282
283
284

        holes = self._holes
        particles = self._particles

        occA = self._occA
Jacob Davison's avatar
Jacob Davison committed
285
        occA2 = self._occA2
Davison's avatar
Davison committed
286
287
        occB = self._occB
        occC = self._occC
288
        occD = self._occD
Davison's avatar
Davison committed
289
290
291
292
293
        occF = self._occF
        occG = self._occG
        occH = self._occH
        occI = self._occI
        occJ = self._occJ
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313

        # Calculate 1B generator
        # fourth term
        sum4_1b_1 = np.matmul(np.transpose(occD,[2,3,0,1]), God)
        sum4_1b_2 = np.matmul(np.transpose(occD,[2,3,0,1]), Gd)
        sum4_1b_3 = ncon([Wd,  sum4_1b_1], [(0,1,-1,2,3,-2),(2,3,0,1)]).numpy()
        sum4_1b_4 = ncon([Wod, sum4_1b_2], [(0,1,-1,2,3,-2),(2,3,0,1)]).numpy()
        sum4_1b = sum4_1b_3 - sum4_1b_4

        # fifth term
        sum5_1b_1 = ncon([occF, Wd],
                         [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
        sum5_1b_2 = ncon([occF, Wod],
                         [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)]).numpy()
        sum5_1b_3 = ncon([sum5_1b_1, Wod],
                         [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
        sum5_1b_4 = ncon([sum5_1b_2, Wd],
                         [(0,1,-1,2,3,4), (2,3,4,0,1,-2)]).numpy()
        sum5_1b = sum5_1b_3 - sum5_1b_4

Jacob Davison's avatar
Jacob Davison committed
314
        eta1B += 0.25*sum4_1b + (1/12)*sum5_1b
315

Davison's avatar
Davison committed
316
317
        # Calculate 2B generator
        # fourth term
Davison's avatar
Davison committed
318
319
        sum4_2b_1 = np.matmul(-1.0*np.transpose(occA2), fod)
        sum4_2b_2 = np.matmul(-1.0*np.transpose(occA2),  fd)
Jacob Davison's avatar
Jacob Davison committed
320
321
        sum4_2b_3 = ncon([Wd,  sum4_2b_1], [(0,-1,-2,1,-3,-4), (1,0)]).numpy()
        sum4_2b_4 = ncon([Wod, sum4_2b_2], [(0,-1,-2,1,-3,-4), (1,0)]).numpy()
Davison's avatar
Davison committed
322
323
324
        sum4_2b = sum4_2b_3 - sum4_2b_4

        #fifth term
Jacob Davison's avatar
Jacob Davison committed
325
326
327
328
329
        sum5_2b_1 = ncon([occG, God], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
        sum5_2b_2 = ncon([occG,  Gd], [(-1,-2,-4,0,1,2), (1,2,0,-3)]).numpy()
        sum5_2b_3 = ncon([Wd,  sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
        sum5_2b_4 = ncon([Wod, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)]).numpy()
        sum5_2b_5 = sum5_2b_3 - sum5_2b_4
Jacob Davison's avatar
Jacob Davison committed
330
        sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
Jacob Davison's avatar
Jacob Davison committed
331
332
333
334
335
336
337
338
                    np.transpose(sum5_2b_5, [0,1,3,2]) + \
                    np.transpose(sum5_2b_5, [2,3,0,1])

        #sixth term
        sum6_2b_1 = ncon([occH, Wod], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)]).numpy()
        sum6_2b_2 = ncon([occH, Wod], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)]).numpy()
        sum6_2b_3 = ncon([Wd, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)]).numpy()
        sum6_2b_4 = ncon([Wd, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)]).numpy()
Jacob Davison's avatar
Jacob Davison committed
339
        sum6_2b = sum6_2b_3 - sum6_2b_4
Jacob Davison's avatar
Jacob Davison committed
340
341

        #seventh term
Davison's avatar
Davison committed
342
        sum7_2b_1 = ncon([occI, Wod], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)]).numpy()
Jacob Davison's avatar
Jacob Davison committed
343
344
345
346
        sum7_2b_2 = ncon([Wd, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)]).numpy()
        sum7_2b = sum7_2b_2 - np.transpose(sum7_2b_2,[1,0,2,3]) - \
                              np.transpose(sum7_2b_2,[0,1,3,2]) + \
                              np.transpose(sum7_2b_2,[1,0,3,2])
Jacob Davison's avatar
Jacob Davison committed
347

Jacob Davison's avatar
Jacob Davison committed
348
        eta2B += sum4_2b + 0.5*sum5_2b + (1/6)*sum6_2b + 0.25*sum7_2b
Davison's avatar
Davison committed
349

Jacob Davison's avatar
Jacob Davison committed
350
        return (eta1B, eta2B)