generator.py 8.59 KB
Newer Older
1
2
from tensornetwork import *
import numpy as np
Davison's avatar
Davison committed
3
4
from oop_imsrg.hamiltonian import *
from oop_imsrg.occupation_tensors import *
5

6
7
8
9
10
11
12
13
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
14
15
    """Calculate Wegner's generator for a normal ordered Hamiltonian.
       Truncate at two-body interactions."""
16
17

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

        Arguments:

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

25
        assert isinstance(h, Hamiltonian), "Arg 0 must be Hamiltonian object"
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
        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):
41
42
43
        """Returns:

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

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

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

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

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

Davison's avatar
Davison committed
63
    def decouple_OD(self):
Davison's avatar
Davison committed
64
65
        """Decouple the off-/diagonal elements from each other in
        the one- and two-body tensors. This procedure is outlined in
66
67
68
69
70
71
72
73
74
        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"""

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
        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)

93
    def calc_eta(self):
Davison's avatar
Davison committed
94
        """Calculate the generator. The terms are defined in An
95
96
97
98
99
100
101
102
        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"""

Davison's avatar
Davison committed
103
        fd, fod, Gd, God = self.decouple_OD()
104
105
106
107
108

        holes = self._holes
        particles = self._particles

        occA = self._occA
109
        occB = self._occB
110
111
112
113
114
115
116
        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
117

118
119
120
121
122
        # 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
123

124
125
126
127
128
        # 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
129

130
131
132
133
134
135
136
137
138
        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
139

140
141
142
143
144
145
        # 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
146

147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
        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

166
167
        return (eta1B, eta2B)

Davison's avatar
Davison committed
168
class WegnerGenerator3B(WegnerGenerator):
169

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

Davison's avatar
Davison committed
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
    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
        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

    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"""
229

Davison's avatar
Davison committed
230
        fd, fod, Gd, God = super().decouple_OD()
Davison's avatar
Davison committed
231

Davison's avatar
Davison committed
232
        W = self.W
233

Davison's avatar
Davison committed
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
        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)

    def calc_eta(self):
        """Inherits from WegnerGenerator.

        Calculate the generator. The terms are defined in An
        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
         eta3B) -- three-body generator"""

Davison's avatar
Davison committed
257
        eta1B, eta2B = super().calc_eta()
Davison's avatar
Davison committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271

        fd, fod, Gd, God, Wd, Wod = self.decouple_OD()

        holes = self._holes
        particles = self._particles

        occA = self._occA
        occB = self._occB
        occC = self._occC
        occF = self._occF
        occG = self._occG
        occH = self._occH
        occI = self._occI
        occJ = self._occJ