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

9
10
from numba import jit

11
12
13
14
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."""
15
    
16
17
18
19
    def calc_eta():
        print("Function that calculates the generator")

class WegnerGenerator(Generator):
Davison's avatar
Davison committed
20
21
    """Calculate Wegner's generator for a normal ordered Hamiltonian.
       Truncate at two-body interactions."""
22
23

    def __init__(self, h, occ_t):
24
25
26
27
28
29
        """Class constructor. Instantiate WegnerGenerator object.

        Arguments:

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

Davison, Jacob's avatar
Davison, Jacob committed
31
#        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
32
33
34
35
36
        assert isinstance(occ_t, OccupationTensors), "Arg 1 must be OccupationTensors object"

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

37

38
39
40
41
        self._holes = h.holes
        self._particles = h.particles

        self._occA = occ_t.occA
42
        self._occA4 = occ_t.occA4
43
        self._occB = occ_t.occB
44
        self._occB4 = occ_t.occB4
45
46
47
        self._occC = occ_t.occC
        self._occD = occ_t.occD

48
49
50
51
52
53
54
        ref = h.reference
        n = len(self._holes)+len(self._particles)
        Ga = tn.Node(np.transpose(np.append(ref[np.newaxis,:], np.zeros((1,n)), axis=0).astype(float)))
        Gb = tn.Node(np.append(ref[::-1][np.newaxis,:],np.zeros((1,n)), axis=0).astype(float))
        
        self._occRef1 = tn.ncon([Ga,Gb], [(-1,1),(1,-2)])                                                                    # n_a(1-n_b)
        self._occRef2 = tn.ncon([tn.Node(np.transpose(Gb.tensor)), tn.Node(np.transpose(Ga.tensor))], [(-1,1),(1,-2)])       # (1-n_a)n_b
55

56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
        self._eta1B = np.zeros_like(self.f)
        self._eta2B = np.zeros_like(self.G)

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

        eta1B -- one-body generator"""
        return self._eta1B

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

        eta2B -- one-body generator"""
        return self._eta2B


74
75
    @property
    def f(self):
76
77
78
        """Returns:

        f -- one-body tensor elements (initialized by Hamiltonian object)"""
79
80
81
82
        return self._f

    @property
    def G(self):
83
84
        """Returns:

Davison's avatar
Davison committed
85
        f -- two-body tensor elements (initialized by Hamiltonian object)"""
86
87
88
89
        return self._G

    @f.setter
    def f(self, f):
90
        """Sets the one-body tensor."""
91
92
93
94
        self._f = f

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

98
    # @classmethod
Davison's avatar
Davison committed
99
    def decouple_OD(self):
Davison's avatar
Davison committed
100
101
        """Decouple the off-/diagonal elements from each other in
        the one- and two-body tensors. This procedure is outlined in
102
103
104
105
106
107
108
109
110
        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"""

111
112
113
114
115
        f = self.f
        G = self.G
        holes = self._holes
        particles = self._particles

116
117
118
119
120
121
        occRef1 = self._occRef1
        occRef2 = self._occRef2

        occD = self._occD
        occDt = tn.Node(np.transpose(occD.tensor,[2,3,0,1]))

122
        # - Decouple off-diagonal 1B and 2B pieces
123
 #       fod1 = tn.ncon([])
124
        fod = np.zeros(f.shape, dtype=np.float32)
125
126
127
128
        fod += np.multiply(occRef2.tensor, f)
        fod += np.multiply(occRef1.tensor, f)
        # fod[np.ix_(particles, holes)] += f[np.ix_(particles, holes)]
        # fod[np.ix_(holes, particles)] += f[np.ix_(holes, particles)]
129
        fd = f - fod
130

131
        God = np.zeros(G.shape, dtype=np.float32)
132
133
134
135
        God += np.multiply(occDt.tensor, G)
        God += np.multiply(occD.tensor, G)
        # 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)]
136
        Gd = G - God
137

138
        return (fd, fod, Gd, God)
139

140
    # @classmethod
141
    def calc_eta(self):
Davison's avatar
Davison committed
142
        """Calculate the generator. The terms are defined in An
143
144
145
146
147
148
149
150
        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"""

151
        partition = self.decouple_OD()
152
153
154
155
        fd = partition[0].astype(np.float32)
        fod = partition[1].astype(np.float32)
        Gd = partition[2].astype(np.float32)
        God = partition[3].astype(np.float32)
156
157
158
159
160

        holes = self._holes
        particles = self._particles

        occA = self._occA
161
        occA4 = self._occA4
162
        occB = self._occB
163
        occB4 = self._occB4
164
165
        occC = self._occC

166

167
168
        # - Calculate 1B generator
        # first term
169
        sum1_1b_1 = tn.ncon([fd, fod], [(-1, 1), (1, -2)])#.numpy()
170
        sum1_1b_2 = np.transpose(sum1_1b_1)
171
        sum1_1b = sum1_1b_1 - sum1_1b_2
Davison's avatar
Davison committed
172

173
        # second term
174
175
176
177
178
179
180
181
        # sum2_1b_1 = tn.ncon([fd, occA, God], [(3,4), (3,4,1,2), (2,-1,1,-2)])#.numpy()
        # sum2_1b_2 = tn.ncon([fod, occA, Gd], [(3,4), (3,4,1,2), (2,-1,1,-2)])#.numpy()
        # sum2_1b = sum2_1b_1 - sum2_1b_2

        fdPrime = np.multiply(occA.tensor, fd)
        fodPrime = np.multiply(occA.tensor, fod)
        sum2_1b_1 = tn.ncon([fdPrime, God], [(1,2), (2,-1,1,-2)])
        sum2_1b_2 = tn.ncon([fodPrime, Gd], [(1,2), (2,-1,1,-2)])
182
        sum2_1b = sum2_1b_1 - sum2_1b_2
183
        
184
185
        # sum2_1b_1 = tn.ncon([fd, God], [(0, 1), (1, -1, 0, -2)])#.numpy()
        # sum2_1b_2 = tn.ncon([fod, Gd], [(0, 1), (1, -1, 0, -2)])#.numpy()
186
        # sum2_1b_3 = sum2_1b_1 - sum2_1b_2
187
        # sum2_1b = tn.ncon([occA, sum2_1b_3],[(-1, -2, 0, 1), (0,1)])#.numpy()
Davison's avatar
Davison committed
188

189
        # third term
190
191
        #sum3_1b_1 = tn.ncon([Gd, occC, God], [(6,-1,4,5), (4,5,6,1,2,3), (1,2,3,-2)])#.numpy()
        #sum3_1b = sum3_1b_1 - np.transpose(sum3_1b_1)
192
        
193
194
195
196
        sum3_1b_1 = np.multiply(occC.tensor, God)#np.multiply(tn.outer_product(tn.Node(occC), tn.Node(np.ones(8))).tensor, God)
        sum3_1b_2 = tn.ncon([Gd, God], [(3,-1,1,2),(1,2,3,-2)])
        sum3_1b = sum3_1b_2 - np.transpose(sum3_1b_2)

197
198
        # sum3_1b_1 = tn.ncon([occC, God], [(-1, -2, -3, 0, 1, 2), (0, 1, 2, -4)])#.numpy()
        # sum3_1b_2 = tn.ncon([Gd, sum3_1b_1], [(2, -1, 0, 1), (0, 1, 2, -2)])#.numpy()
199
200
        # sum3_1b_3 = np.transpose(sum3_1b_2)
        # sum3_1b = sum3_1b_2 - sum3_1b_3
Davison's avatar
Davison committed
201

202
        eta1B = sum1_1b + sum2_1b + 0.5*sum3_1b
203
204
205

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

212
        # first term (P_kl piece)
213
214
        sum1_2b_6 = tn.ncon([fd, God], [(1, -3), (-1, -2, 1, -4)])#.numpy()
        sum1_2b_7 = tn.ncon([fod, Gd], [(1, -3), (-1, -2, 1, -4)])#.numpy()
215
        sum1_2b_8 = sum1_2b_6 - sum1_2b_7
216
        sum1_2b_9 = np.transpose(sum1_2b_8, [0, 1, 3, 2])
217
        sum1_2b_10 = sum1_2b_8 - sum1_2b_9
Davison's avatar
Davison committed
218

219
220
221
        sum1_2b = sum1_2b_5 - sum1_2b_10

        # second term
222
223
        #sum2_2b_1 = tn.ncon([Gd, occB, God], [(-1,-2,3,4), (3,4,1,2), (1,2,-3,-4)])#.numpy()
        #sum2_2b_2 = tn.ncon([God, occB, Gd], [(-1,-2,3,4), (3,4,1,2), (1,2,-3,-4)])#.numpy()
224

225
226
227
228
229
        
        GodPrime = np.multiply(occB4.tensor, God)
        GdPrime = np.multiply(occB4.tensor, Gd)
        sum2_2b_1 = tn.ncon([Gd, GodPrime], [(-1,-2,1,2), (1,2,-3,-4)])
        sum2_2b_2 = tn.ncon([God, GdPrime], [(-1,-2,1,2), (1,2,-3,-4)])
230
231
232
233

#        sum2_2b_1 = tn.ncon([Gd, occB[:,:,None],occB[None,:,:], God], [(-1,-2,4,5), (4,5,1), (1,2,3), (2,3,-3,-4)])#.numpy()
#        sum2_2b_2 = tn.ncon([God, occB[:,:,None],occB[None,:,:], Gd], [(-1,-2,4,5), (4,5,1), (1,2,3), (2,3,-3,-4)])#.numpy()

234
        sum2_2b = sum2_2b_1 - sum2_2b_2
235

236
237
238
239
        # sum2_2b_1 = tn.ncon([occB, God], [(-1, -2, 0, 1), (0, 1, -3, -4)])#.numpy()
        # sum2_2b_2 = tn.ncon([occB,  Gd], [(-1, -2, 0, 1), (0, 1, -3, -4)])#.numpy()
        # sum2_2b_3 = tn.ncon([Gd,  sum2_2b_1], [(-1, -2, 0, 1), (0, 1, -3, -4)])#.numpy()
        # sum2_2b_4 = tn.ncon([God, sum2_2b_2], [(-1, -2, 0, 1), (0, 1, -3, -4)])#.numpy()
240
        # sum2_2b = sum2_2b_3 - sum2_2b_4
241
242

        # third term
243
244
245
246
247
248
        # sum3_2b_1 = tn.ncon([Gd, occA, God], [(3,-1,4,-3), (3,4,1,2), (2,-2,1,-4)])#.numpy()
        # sum3_2b_2 = sum3_2b_1 - np.transpose(sum3_2b_1, [0,1,3,2])
        # sum3_2b = sum3_2b_2 - np.transpose(sum3_2b_2, [1,0,2,3])

        GodPrime = np.multiply(np.transpose(occA4.tensor, [0,2,1,3]), God)
        sum3_2b_1 = tn.ncon([Gd, GodPrime], [(2,-2,1,-4), (1,-1,2,-3)])
249
250
        sum3_2b_2 = sum3_2b_1 - np.transpose(sum3_2b_1, [0,1,3,2])
        sum3_2b = sum3_2b_2 - np.transpose(sum3_2b_2, [1,0,2,3])
251

252
        # sum3_2b_1 = tn.ncon([Gd, God], [(0, -1, 1, -3), (1, -2, 0, -4)])#.numpy()
253
254
255
256
        # 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
257
        # sum3_2b = tn.ncon([occA, sum3_2b_5], [(0, 1, -1, -2), (0, 1, -3, -4)])#.numpy()
258

259
        eta2B = sum1_2b + 0.5*sum2_2b - sum3_2b
260

261
262
        return (eta1B, eta2B)

Davison's avatar
Davison committed
263
264
265
class WegnerGenerator3B(WegnerGenerator):
    """Calculate Wegner's generator for a normal ordered Hamiltonian.
       Truncate at three-body interactions. Inherits from WegnerGenerator."""
266

Davison's avatar
Davison committed
267
268
269
270
271
272
273
274
275
276
277
278
279
    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
280
        self.W = np.zeros(h.n_sp_states*np.ones(6,dtype=np.int32),dtype=np.float32)
Davison's avatar
Davison committed
281
282
283
284
285

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

        self._occA = occ_t.occA
Jacob Davison's avatar
Jacob Davison committed
286
        self._occA2 = occ_t.occA2
Davison's avatar
Davison committed
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
        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

309
    # @classmethod
Davison's avatar
Davison committed
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
    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"""
325

326
327
328
329
330
        partition = super().decouple_OD()
        fd = partition[0]
        fod = partition[1]
        Gd = partition[2]
        God = partition[3]
331
        
Davison's avatar
Davison committed
332
        W = self.W
333
        
Davison's avatar
Davison committed
334
335
336
        holes = self._holes
        particles = self._particles

337
        Wod = np.zeros(W.shape, dtype=np.float32)
Davison's avatar
Davison committed
338
339
340
341
342
343
        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)

344
    # @classmethod
Davison's avatar
Davison committed
345
346
347
    def calc_eta(self):
        """Inherits from WegnerGenerator.

348
349
        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
350
351
352
353
354
355
356

        Returns:

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

357
358
359
360
361
362
363
        partition = self.decouple_OD()
        fd = partition[0]
        fod = partition[1]
        Gd = partition[2]
        God = partition[3]
        Wd = partition[4]
        Wod = partition[5]
364
        
365
        eta1B, eta2B = super().calc_eta()
Davison's avatar
Davison committed
366
367
368
369
370

        holes = self._holes
        particles = self._particles

        occA = self._occA
Jacob Davison's avatar
Jacob Davison committed
371
        occA2 = self._occA2
Davison's avatar
Davison committed
372
373
        occB = self._occB
        occC = self._occC
374
        occD = self._occD
Davison's avatar
Davison committed
375
376
377
378
379
        occF = self._occF
        occG = self._occG
        occH = self._occH
        occI = self._occI
        occJ = self._occJ
380

381

382
383
        # Calculate 1B generator
        # fourth term
384
385
        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)
386
387
        sum4_1b_3 = tn.ncon([Wd,  sum4_1b_1], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
        sum4_1b_4 = tn.ncon([Wod, sum4_1b_2], [(1,2,-1,3,4,-2),(3,4,1,2)])#.numpy()
388
389
390
        sum4_1b = sum4_1b_3 - sum4_1b_4

        # fifth term
391
392
393
394
395
396
        sum5_1b_1 = tn.ncon([Wd, occF, Wod], [(6,7,-1,8,9,10),
                                              (6,7,8,9,10,1,2,3,4,5),
                                              (3,4,5,1,2,-2)])#.numpy()
        sum5_1b_2 = tn.ncon([Wod, occF, Wd], [(6,7,-1,8,9,10),
                                              (6,7,8,9,10,1,2,3,4,5),
                                              (3,4,5,1,2,-2)])#.numpy()
397
        sum5_1b = sum5_1b_1 - sum5_1b_2
398

399
        # sum5_1b_1 = tn.ncon([occF, Wd.astype(np.float32)],
400
        #                  [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)])#.numpy()
401
        # sum5_1b_2 = tn.ncon([occF, Wod.astype(np.float32)],
402
        #                  [(-1,-3,-4,-5,-6,0,1,2,3,4), (0,1,-2,2,3,4)])#.numpy()
403
        # sum5_1b_3 = tn.ncon([sum5_1b_1, Wod.astype(np.float32)],
404
        #                  [(0,1,-1,2,3,4), (2,3,4,0,1,-2)])#.numpy()
405
        # sum5_1b_4 = tn.ncon([sum5_1b_2, Wd.astype(np.float32)],
406
        #                  [(0,1,-1,2,3,4), (2,3,4,0,1,-2)])#.numpy()
407
        # sum5_1b = sum5_1b_3 - sum5_1b_4
408

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

Davison's avatar
Davison committed
411
412
        # Calculate 2B generator
        # fourth term
413
414
        sum4_2b_1 = np.matmul(-1.0*np.transpose(occA2), fod)
        sum4_2b_2 = np.matmul(-1.0*np.transpose(occA2),  fd)
415
416
        sum4_2b_3 = tn.ncon([Wd,  sum4_2b_1], [(1,-1,-2,2,-3,-4), (2,1)])#.numpy()
        sum4_2b_4 = tn.ncon([Wod, sum4_2b_2], [(1,-1,-2,2,-3,-4), (2,1)])#.numpy()
Davison's avatar
Davison committed
417
418
419
        sum4_2b = sum4_2b_3 - sum4_2b_4

        #fifth term
420
421
422
423
424
425
        sum5_2b_1 = tn.ncon([Wd, occG, God], [(4,-1,-2,5,6,-4),
                                              (4,5,6,1,2,3),
                                              (2,3,1,-3)])#.numpy()
        sum5_2b_2 = tn.ncon([Wod, occG, Gd], [(4,-1,-2,5,6,-4),
                                              (4,5,6,1,2,3),
                                              (2,3,1,-3)])#.numpy()
426
427
428
        sum5_2b = sum5_2b_2 - np.transpose(sum5_2b_2, [3,2,0,1]) - \
                    np.transpose(sum5_2b_2, [0,1,3,2]) + \
                    np.transpose(sum5_2b_2, [2,3,0,1])
429
430
431
432
433

        # sum5_2b_1 = tn.ncon([occG, God], [(-1,-2,-4,0,1,2), (1,2,0,-3)])#.numpy()
        # sum5_2b_2 = tn.ncon([occG,  Gd], [(-1,-2,-4,0,1,2), (1,2,0,-3)])#.numpy()
        # sum5_2b_3 = tn.ncon([Wd,  sum5_2b_1], [(0,-1,-2,1,2,-4), (1,2,0,-3)])#.numpy()
        # sum5_2b_4 = tn.ncon([Wod, sum5_2b_2], [(0,-1,-2,1,2,-4), (1,2,0,-3)])#.numpy()
434
435
436
437
        # sum5_2b_5 = sum5_2b_3 - sum5_2b_4
        # sum5_2b = sum5_2b_5 - np.transpose(sum5_2b_5, [3,2,0,1]) - \
        #             np.transpose(sum5_2b_5, [0,1,3,2]) + \
        #             np.transpose(sum5_2b_5, [2,3,0,1])
Jacob Davison's avatar
Jacob Davison committed
438
439

        #sixth term
440
441
442
443
444
445
        sum6_2b_1 = tn.ncon([Wd, occH, Wod], [(5,-1,-2,6,7,8),
                                              (5,6,7,8,1,2,3,4),
                                              (2,3,4,1,-3,-4)])#.numpy()
        sum6_2b_2 = tn.ncon([Wd, occH, Wod], [(6,7,8,5,-3,-4),
                                              (5,6,7,8,1,2,3,4),
                                              (1,-1,-2,2,3,4)])#.numpy()
446
447
        sum6_2b = sum6_2b_1 - sum6_2b_2

448
449
450
451
        # sum6_2b_1 = tn.ncon([occH, Wod], [(-1,-2,-3,-4,0,1,2,3),(1,2,3,0,-5,-6)])#.numpy()
        # sum6_2b_2 = tn.ncon([occH, Wod], [(-3,-4,-5,-6,0,1,2,3),(0,-1,-2,1,2,3)])#.numpy()
        # sum6_2b_3 = tn.ncon([Wd, sum6_2b_1], [(0,-1,-2,1,2,3), (1,2,3,0,-3,-4)])#.numpy()
        # sum6_2b_4 = tn.ncon([Wd, sum6_2b_2], [(1,2,3,0,-3,-4), (0,-1,-2,1,2,3)])#.numpy()
452
        # sum6_2b = sum6_2b_3 - sum6_2b_4
Jacob Davison's avatar
Jacob Davison committed
453
454

        #seventh term
455
456
457
        sum7_2b_1 = tn.ncon([Wd, occI, Wod], [(5,6,-1,7,8,-4),
                                              (5,6,7,8,1,2,3,4),
                                              (3,4,-2,1,2,-3)])#.numpy()
458
459
460
        sum7_2b = sum7_2b_1 - np.transpose(sum7_2b_1,[1,0,2,3]) - \
                              np.transpose(sum7_2b_1,[0,1,3,2]) + \
                              np.transpose(sum7_2b_1,[1,0,3,2])
461
462
463

        # sum7_2b_1 = tn.ncon([occI, Wod], [(-1,-2,-3,-4,0,1,2,3), (2,3,-5,0,1,-6)])#.numpy()
        # sum7_2b_2 = tn.ncon([Wd, sum7_2b_1], [(0,1,-1,2,3,-4),(2,3,-2,0,1,-3)])#.numpy()
464
465
466
        # 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
467

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

Davison's avatar
Davison committed
470
471
472
473
        # Calculate 3B generator
        #first, second, and third terms (one index contraction)

        #terms with P(i/jk) -- line 1 and 2
474
475
        sum1_3b_1 = tn.ncon([fd, Wod], [(-1,1), (1,-2,-3,-4,-5,-6)])#.numpy()
        sum1_3b_2 = tn.ncon([fod, Wd], [(-1,1), (1,-2,-3,-4,-5,-6)])#.numpy()
Davison's avatar
Davison committed
476
        sum1_3b_3 = sum1_3b_1 - sum1_3b_2
477
478
        sum1_3b_4 = sum1_3b_3 - np.transpose(sum1_3b_3, [1,0,2,3,4,5]) - \
                                np.transpose(sum1_3b_3, [2,1,0,3,4,5])
Davison's avatar
Davison committed
479
480

        #terms with P(l/mn) -- line 1 and 2
481
482
        sum1_3b_5 = tn.ncon([fd, Wod], [(1,-4), (-1,-2,-3,1,-5,-6)])#.numpy()
        sum1_3b_6 = tn.ncon([fod, Wd], [(1,-4), (-1,-2,-3,1,-5,-6)])#.numpy()
Davison's avatar
Davison committed
483
        sum1_3b_7 = sum1_3b_6 - sum1_3b_5
484
485
        sum1_3b_8 = sum1_3b_7 - np.transpose(sum1_3b_7, [0,1,2,4,3,5]) - \
                                np.transpose(sum1_3b_7, [0,1,2,5,4,3])
Davison's avatar
Davison committed
486
487

        #terms with P(ij/k)P(l/mn) -- line 3
488
489
        sum1_3b_9  = tn.ncon([Gd, God], [(-1,-2,-4,1),(1,-3,-5,-6)])#.numpy()
        sum1_3b_10 = tn.ncon([God, Gd], [(-1,-2,-4,1),(1,-3,-5,-6)])#.numpy()
Davison's avatar
Davison committed
490
        sum1_3b_11 = sum1_3b_9 - sum1_3b_10
491
492
493
494
        sum1_3b_12 = sum1_3b_11 - np.transpose(sum1_3b_11, [0,1,2,4,3,5]) - \
                                  np.transpose(sum1_3b_11, [0,1,2,5,4,3])
        sum1_3b_13 = sum1_3b_12 - np.transpose(sum1_3b_12, [2,1,0,3,4,5]) - \
                                  np.transpose(sum1_3b_12, [0,2,1,3,4,5])
Davison's avatar
Davison committed
495
496
497
498

        sum1_3b = sum1_3b_4 + sum1_3b_8 + sum1_3b_13

        #fourth term
499
500
        sum4_3b_1 = tn.ncon([Gd, occB, Wod], [(-1,-2,3,4),(3,4,1,2),(1,2,-3,-4,-5,-6)])#.numpy()
        sum4_3b_2 = tn.ncon([God, occB, Wd], [(-1,-2,3,4),(3,4,1,2),(1,2,-3,-4,-5,-6)])#.numpy()
Davison's avatar
Davison committed
501
        sum4_3b_3 = sum4_3b_1 - sum4_3b_2
502
503
        sum4_3b = sum4_3b_3 - np.transpose(sum4_3b_3, [1,0,2,3,4,5]) - \
                              np.transpose(sum4_3b_3, [2,1,0,3,4,5])
Davison's avatar
Davison committed
504
505

        #fifth term
506
507
        sum5_3b_1 = tn.ncon([Gd, occB, Wod], [(3,4,-4,-5),(3,4,1,2),(-1,-2,-3,1,2,-6)])#.numpy()
        sum5_3b_2 = tn.ncon([God, occB, Wd], [(3,4,-4,-5),(3,4,1,2),(-1,-2,-3,1,2,-6)])#.numpy()
Davison's avatar
Davison committed
508
        sum5_3b_3 = sum5_3b_1 - sum5_3b_2
509
510
        sum5_3b = sum5_3b_3 - np.transpose(sum5_3b_3, [0,1,2,5,4,3]) - \
                              np.transpose(sum5_3b_3, [0,1,2,3,5,4])
Davison's avatar
Davison committed
511
512

        #sixth term
513
514
        sum6_3b_1 = tn.ncon([Gd, occA, Wod], [(4,-1,3,-4),(3,4,1,2),(1,-2,-3,2,-5,-6)])#.numpy()
        sum6_3b_2 = tn.ncon([God, occA, Wd], [(4,-1,3,-4),(3,4,1,2),(1,-2,-3,2,-5,-6)])#.numpy()
Davison's avatar
Davison committed
515
        sum6_3b_3 = sum6_3b_1 - sum6_3b_2
516
517
518
519
        sum6_3b_4 = sum6_3b_3 - np.transpose(sum6_3b_3, [0,1,2,4,3,5]) - \
                                np.transpose(sum6_3b_3, [0,1,2,5,4,3])
        sum6_3b = sum6_3b_4 - np.transpose(sum6_3b_4, [1,0,2,3,4,5]) - \
                              np.transpose(sum6_3b_4, [2,1,0,3,4,5])
Davison's avatar
Davison committed
520
521

        #seventh term
522
523
        sum7_3b_1 = tn.ncon([Wd, occJ, Wod], [(-1,-2,-3,4,5,6), (4,5,6,1,2,3), (1,2,3,-4,-5,-6)])#.numpy()
        sum7_3b_2 = tn.ncon([Wod, occJ, Wd], [(-1,-2,-3,4,5,6), (4,5,6,1,2,3), (1,2,3,-4,-5,-6)])#.numpy()
Davison's avatar
Davison committed
524
525
526
        sum7_3b = sum7_3b_1 - sum7_3b_2

        #eighth term
527
528
        sum8_3b_1 = tn.ncon([Wd, occC, Wod], [(4,5,-3,6,-5,-6), (4,5,6,1,2,3), (3,-1,-2,1,2,-4)])#.numpy()
        sum8_3b_2 = tn.ncon([Wd, occC, Wod], [(6,-2,-3,4,5,-6), (4,5,6,1,2,3), (-1,1,2,-4,-5,3)])#.numpy()
Davison's avatar
Davison committed
529
        sum8_3b_3 = sum8_3b_1 - sum8_3b_2
530
531
532
533
        sum8_3b_4 = sum8_3b_3 - np.transpose(sum8_3b_3, [0,1,2,4,3,5]) - \
                                np.transpose(sum8_3b_3, [0,1,2,5,4,3])
        sum8_3b = sum8_3b_4 - np.transpose(sum8_3b_4, [2,1,0,3,4,5]) - \
                              np.transpose(sum8_3b_4, [0,2,1,3,4,5])
Davison's avatar
Davison committed
534
535
536
537

        eta3B = sum1_3b + 0.5*sum4_3b + (-0.5)*sum5_3b + (-1.0)*sum6_3b + (1/6)*sum7_3b + 0.5*sum8_3b

        return (eta1B, eta2B, eta3B)
538
539
540
541
542
543


class WhiteGenerator(Generator):
    """Calculate White's generator for a normal ordered Hamiltonian.
       This standard implemenation uses Epstein-Nesbet denominators."""

544
    
545
546
    def __init__(self, h):

Davison, Jacob's avatar
Davison, Jacob committed
547
#        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
548
        
549
550
551
        self._my_f = h.f
        self._my_G = h.G
        
552
553
554
555
        self._holes = h.holes
        self._particles = h.particles
        self._sp_basis = h.sp_basis

556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
        self._eta1B = np.zeros_like(self.f)
        self._eta2B = np.zeros_like(self.G)

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

        eta1B -- one-body generator"""
        return self._eta1B

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

        eta2B -- one-body generator"""
        return self._eta2B


574
575
576
577
578
    @property
    def f(self):
        """Returns:

        f -- one-body tensor elements (initialized by Hamiltonian object)"""
579
        return self._my_f
580
581
582
583
584
585

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

        f -- two-body tensor elements (initialized by Hamiltonian object)"""
586
        return self._my_G
587
588
589
590

    @f.setter
    def f(self, f):
        """Sets the one-body tensor."""
591
        self._my_f = f
592
593
594
595

    @G.setter
    def G(self, G):
        """Sets the two-body tensor."""
596
597
        self._my_G = G

598
599
600
601
602
603
604
605
606

    def calc_eta(self):

        bas1B = self._sp_basis
        holes = self._holes
        particles = self._particles

        f = self.f
        G = self.G
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623

        eta1B, eta2B = self._wrapper_calc_eta(bas1B, holes, particles, f, G)
        self._eta1B = eta1B
        self._eta2B = eta2B

        return (eta1B, eta2B)

    @staticmethod
    @jit(nopython=True)
    def _wrapper_calc_eta(bas1B, holes, particles, f, G):

        # bas1B = self._sp_basis
        # holes = self._holes
        # particles = self._particles

        # f = self.f
        # G = self.G
624
625
626
627
628
629
630
        
        eta1B = np.zeros_like(f)
        eta2B = np.zeros_like(G)
        
        for a in particles:
            for i in holes:
                denom = f[a,a] - f[i,i] + G[a,i,a,i]
631
632
633
634
635
636

                if abs(denom)<1.0e-10:
                    result = 0.25 * np.pi * np.sign(f[a,i]) * np.sign(denom)
                else:
                    result = f[a,i] / denom

637
638
639
640

                eta1B[a,i] = result
                eta1B[i,a] = -result
                
641
642
                # if denom < 1:
                #     print('one body {}{},'.format(a, i), denom)
643
644
645
646
647
648
649
650
651
652
653
654
                # if a == 4 and i == 3:
                #     print(G[a,i,a,i])

        for a in particles:
            for b in particles:
                for i in holes:
                    for j in holes:
                        denom = (
                            f[a,a] + f[b,b] - f[i,i] - f[j,j] 
                            + G[a,b,a,b] + G[i,j,i,j] - G[a,i,a,i]
                            - G[b,j,b,j] - G[a,j,a,j] - G[b,i,b,i]
                        )
655
656
657
658
659
660
661

                        if abs(denom)<1.0e-10:
                            result = 0.25 * np.pi * np.sign(G[a,b,i,j]) * np.sign(denom)
                        else:
                            result = G[a,b,i,j] / denom


662
663
664
665
666
667
668
669
670
                        
                        eta2B[a,b,i,j] = result
                        eta2B[i,j,a,b] = -result

                        # if denom < 1:
                        #    print('two body {}{}{}{},'.format(a,b,i,j), denom)
                        # if a == 5 and b == 4 and i == 3 and j ==2:
                        #     print(denom)
                        #print(eta2B[5,4,3,2])
671
672
                    
                    
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
        # # Obtain denominator terms
        # fpp = f[np.ix_(particles), np.ix_(particles)]
        # fhh = f[np.ix_(holes), np.ix_(holes)]
        # Gphph = G[np.ix_(particles), np.ix_(holes), np.ix_(particles), np.ix_(holes)]
        # Gpppp = G[np.ix_(particles), np.ix_(particles), np.ix_(particles), np.ix_(particles)]
        # Ghhhh = G[np.ix_(holes), np.ix_(holes), np.ix_(holes), np.ix_(holes)]

        # # Compute one body denominators (Epstein-Nesbet)
        # denom1B = np.ones_like(f)
        # denom1B[np.ix_(particles), np.ix_(holes)] = fpp - fhh + Gphph
        # #denom1B[np.ix_(holes), np.ix_(particles)] = -denom1B[np.ix_(particles), np.ix_(holes)]
        
        # # Compute two body denominators (Epstein-Nesbet)
        # denom2B = np.ones_like(G)
        # denom2B[np.ix_(particles), np.ix_(particles), np.ix_(holes), np.ix_(holes)] = fpp + fpp - fhh - fhh - Gpppp + Ghhhh - Gphph - Gphph - Gphph - Gphph
        # #denom2B[np.ix_(holes), np.ix_(holes), np.ix_(particles), np.ix_(particles)] = -denom2B[np.ix_(particles), np.ix_(particles), np.ix_(holes), np.ix_(holes)]
        
        # # Compute generator
        # eta1B = np.zeros_like(f)
        # eta1B[np.ix_(particles), np.ix_(holes)] = np.divide(f[np.ix_(particles), np.ix_(holes)], denom1B[np.ix_(particles), np.ix_(holes)])
        # eta1B[np.ix_(holes), np.ix_(particles)] = -eta1B[np.ix_(particles), np.ix_(holes)]
        
        # eta2B = np.zeros_like(G)
        # eta2B[np.ix_(particles), np.ix_(particles), np.ix_(holes), np.ix_(holes)] = np.divide(G[np.ix_(particles), np.ix_(particles), np.ix_(holes), np.ix_(holes)], 
        #                                                                                       denom2B[np.ix_(particles), np.ix_(particles), np.ix_(holes), np.ix_(holes)])
        # eta2B[np.ix_(holes), np.ix_(holes), np.ix_(particles), np.ix_(particles)] = -eta2B[np.ix_(particles), np.ix_(particles), np.ix_(holes), np.ix_(holes)]

        return (eta1B, eta2B)
        
class WhiteGeneratorMP(Generator):
    """Calculate White's generator for a normal ordered Hamiltonian.
       This "standard" implemenation uses Moller-Plesset denominators."""


    def __init__(self, h):

        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
        
        self.f = h.f
        self.G = h.G

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

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

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

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

        f -- two-body tensor elements (initialized by Hamiltonian object)"""
        return self._G

    @f.setter
    def f(self, f):
        """Sets the one-body tensor."""
        self._f = f

    @G.setter
    def G(self, G):
        """Sets the two-body tensor."""
        self._G = G

    def calc_eta(self):

        bas1B = self._sp_basis
        holes = self._holes
        particles = self._particles

        f = self.f
        G = self.G
        
        eta1B = np.zeros_like(f)
        eta2B = np.zeros_like(G)
        
        for a in particles:
            for i in holes:
                denom = f[a,a] - f[i,i]
                result = f[a,i] / denom

                eta1B[a,i] = result
                eta1B[i,a] = -result
                
Jacob August Davison's avatar
Jacob August Davison committed
762
763
                # if denom < 1:
                #     print('one body {}{},'.format(a, i), denom)
764
765
766
767
768
769
770
771
772
773
774
775
776

        for a in particles:
            for b in particles:
                for i in holes:
                    for j in holes:
                        denom = (
                            f[a,a] + f[b,b] - f[i,i] - f[j,j] 
                        )
                        result = G[a,b,i,j] / denom
                        
                        eta2B[a,b,i,j] = result
                        eta2B[i,j,a,b] = -result

Jacob August Davison's avatar
Jacob August Davison committed
777
778
                        # if denom < 1:
                        #     print('two body {}{}{}{},'.format(a,b,i,j), denom)
779
780
781
                        #print(eta2B[5,4,3,2])


782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
        return (eta1B, eta2B)

class BrillouinGenerator(Generator):
    """Calculate Brillouin generator for a normal ordered Hamiltonian."""


    def __init__(self, h):

        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
        
        self.f = h.f
        self.G = h.G

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

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

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

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

        f -- two-body tensor elements (initialized by Hamiltonian object)"""
        return self._G

    @f.setter
    def f(self, f):
        """Sets the one-body tensor."""
        self._f = f

    @G.setter
    def G(self, G):
        """Sets the two-body tensor."""
        self._G = G

    def calc_eta(self):

        bas1B = self._sp_basis
        holes = self._holes
        particles = self._particles

        f = self.f
        G = self.G
        
        eta1B = np.zeros_like(f)
        eta2B = np.zeros_like(G)
        

        for a in particles:
            for i in holes:
                # (1-n_a)n_i - n_a(1-n_i) = n_i - n_a
                eta1B[a, i] =  f[a,i]
                eta1B[i, a] = -f[a,i]



        for a in particles:
            for b in particles:
                for i in holes:
                    for j in holes:
                        val = G[a,b,i,j]

                        eta2B[a,b,i,j] = val
                        eta2B[i,j,a,b] = -val


854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
        return (eta1B, eta2B)

class ImTimeGenerator(Generator):
    """Calculate Imaginary time generator for a normal ordered Hamiltonian."""


    def __init__(self, h):

        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
        
        self.f = h.f
        self.G = h.G

        self._holes = h.holes
        self._particles = h.particles
        self._sp_basis = h.sp_basis
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
        self._eta1B = np.zeros_like(self.f)
        self._eta2B = np.zeros_like(self.G)

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

        eta1B -- one-body generator"""
        return self._eta1B

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

        eta2B -- one-body generator"""
        return self._eta2B
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951

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

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

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

        f -- two-body tensor elements (initialized by Hamiltonian object)"""
        return self._G

    @f.setter
    def f(self, f):
        """Sets the one-body tensor."""
        self._f = f

    @G.setter
    def G(self, G):
        """Sets the two-body tensor."""
        self._G = G

    def calc_eta(self):

        bas1B = self._sp_basis
        holes = self._holes
        particles = self._particles

        f = self.f
        G = self.G
        
        eta1B = np.zeros_like(f)
        eta2B = np.zeros_like(G)
        

        for a in particles:
            for i in holes:
                dE = f[a,a] - f[i,i] + G[a,i,a,i]
                val = np.sign(dE)*f[a,i]
                eta1B[a, i] =  val
                eta1B[i, a] = -val 


        for a in particles:
            for b in particles:
                for i in holes:
                    for j in holes:
                        dE = ( 
                            f[a,a] + f[b,b] - f[i,i] - f[j,j]  
                            + G[a,b,a,b] 
                            + G[i,j,i,j]
                            - G[a,i,a,i] 
                            - G[a,j,a,j] 
                            - G[b,i,b,i] 
                            - G[b,j,b,j] 
                        )

                        val = np.sign(dE)*G[a,b,i,j]

                        eta2B[a,b,i,j] = val
                        eta2B[i,j,a,b] = -val


952
        return (eta1B, eta2B)