main.py 6.96 KB
Newer Older
Davison's avatar
Davison committed
1
# Main program for IM-SRG(2).
2
3
4
5
6
7
8
9
10
11
12


# Author: Jacob Davison
# Date:   07/10/2019

# import packages, libraries, and modules
# libraries
from scipy.integrate import odeint, ode
import numpy as np
import time
import pickle
13
#import tensorflow as tf
Jacob August Davison's avatar
Jacob August Davison committed
14
#tf.enable_v2_behavior()
15
#print("GPU available: ",tf.test.is_gpu_available())
16
17
import tracemalloc
import os, sys
Jacob August Davison's avatar
Jacob August Davison committed
18
#from memory_profiler import profile
19
20
import itertools
import random
21
22
23
24
#import tensornetwork as tn
#tn.set_default_backend("tensorflow")
#sess = tf.compat.v1.Session(config=tf.compat.v1.ConfigProto(log_device_placement=True))
#sess.close()
25

26
27
28
29
30
31
32
# user files
# sys.path.append('C:\\Users\\davison\\Research\\exact_diagonalization\\')
from oop_imsrg.hamiltonian import *
from oop_imsrg.occupation_tensors import *
from oop_imsrg.generator import *
from oop_imsrg.flow import *
from oop_imsrg.plot_data import *
Jacob Davison's avatar
Jacob Davison committed
33
# from oop_imsrg.display_memory import *
34
import oop_imsrg.ci_pairing.cipy_pairing_plus_ph as ci_matrix
35
from oop_imsrg.tests2B import *
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

# @profile
def derivative(t, y, hamiltonian, occ_tensors, generator, flow):
    """Defines the derivative to pass into ode object.

    Arguments:
    (required by scipy.integrate.ode)
    t -- points at which to solve for y
    y -- in this case, 1D array that contains E, f, G

    (additional parameters)
    hamiltonian -- Hamiltonian object
    occ_tensors -- OccupationTensors object
    generator -- Generator object
    flow -- Flow object

    Returns:

    dy -- next step in flow"""

    assert isinstance(hamiltonian, Hamiltonian), "Arg 2 must be Hamiltonian object"
    assert isinstance(occ_tensors, OccupationTensors), "Arg 3 must be OccupationTensors object"
    assert isinstance(generator, Generator), "Arg 4 must be Generator object"
    assert isinstance(flow, Flow), "Arg 5 must be a Flow object"

    E, f, G = ravel(y, hamiltonian.n_sp_states)

63
#    timefi = time.time()
64
    generator.f = f
65
#    timeff = time.time()
66

67
#    timegi = time.time()
68
    generator.G = G
69
70
71
#    timegf = time.time()

#    timefli = time.time()
72
    dE, df, dG = flow.flow(generator)
73
74
75
76
77
78
79
#    timeflf = time.time()

    # print("""
    # f time:     {:2.4f} s
    # G time:     {:2.4f} s
    # Flow time:  {:2.4f} s
    # """.format(timeff-timefi, timegf-timegi, timeflf-timefli))
80

81
    dy = unravel(dE, df, dG)
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

    return dy

# @profile
def unravel(E, f, G):
    """Transforms E, f, and G into a 1D array. Facilitates
    compatability with scipy.integrate.ode.

    Arguments:

    E, f, G -- normal-ordered pieces of Hamiltonian

    Returns:

    concatenation of tensors peeled into 1D arrays"""
    unravel_E = np.reshape(E, -1)
    unravel_f = np.reshape(f, -1)
    unravel_G = np.reshape(G, -1)

    return np.concatenate([unravel_E, unravel_f, unravel_G], axis=0)

# @profile
def ravel(y, bas_len):
    """Transforms 1D array into E, f, and G. Facilitates
    compatability with scipy.integrate.ode.

    Arugments:

    y -- 1D data array (output from unravel)
    bas_len -- length of single particle basis

    Returns:

    E, f, G -- normal-ordered pieces of Hamiltonian"""

    # bas_len = len(np.append(holes,particles))

    ravel_E = np.reshape(y[0], ())
    ravel_f = np.reshape(y[1:bas_len**2+1], (bas_len, bas_len))
    ravel_G = np.reshape(y[bas_len**2+1:bas_len**2+1+bas_len**4],
                         (bas_len, bas_len, bas_len, bas_len))
123
    
124
125
126
127

    return(ravel_E, ravel_f, ravel_G)

# @profile
128
def main(n_holes, n_particles, ref=None, d=1.0, g=0.5, pb=0.0, verbose=1):
Davison's avatar
Davison committed
129
    """Main method uses scipy.integrate.ode to solve the IMSRG(2) flow
130
131
    equations."""

Jacob August Davison's avatar
Jacob August Davison committed
132
133
134
135
    start = time.time() # start full timer

    initi = time.time() # start instantiation timer

136
137
    if ref == None:
        ha = PairingHamiltonian2B(n_holes, n_particles, d=d, g=g, pb=pb)
138
        ref = ha.reference # this is just for printing
139
140
    else:
        ha = PairingHamiltonian2B(n_holes, n_particles, ref=ref, d=d, g=g, pb=pb)
Jacob August Davison's avatar
Jacob August Davison committed
141

142
143
    ot = OccupationTensors(ha.sp_basis, ha.reference)
    wg = WegnerGenerator(ha, ot)
Jacob August Davison's avatar
Jacob August Davison committed
144
145
146
147
    fl = Flow_IMSRG2(ha, ot) 

    initf = time.time() # finish instantiation timer

148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    if verbose:
        print("Initialized objects in {:2.4f} seconds\n".format(initf-initi))

    if verbose:
        print("""Pairing model IM-SRG(2) flow:
        d              = {:2.4f}
        g              = {:2.4f}
        pb             = {:2.4f}
        SP basis size  = {:2d}
        n_holes        = {:2d}
        n_particles    = {:2d}
        ref            = {d}""".format(ha.d, ha.g, ha.pb, ha.n_sp_states,
                                       len(ha.holes), len(ha.particles),
                                       d=ref) )
    if verbose:
        print("Flowing...")
Jacob August Davison's avatar
Jacob August Davison committed
164
    flowi = time.time()
165
166

    # --- Solve the IM-SRG flow
167
168
    y0 = unravel(ha.E, ha.f, ha.G)

169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
    solver = ode(derivative,jac=None)
    solver.set_integrator('vode', method='bdf', order=5, nsteps=500)
    solver.set_f_params(ha, ot, wg, fl)
    solver.set_initial_value(y0, 0.)

    sfinal = 50
    ds = 0.1
    s_vals = []
    E_vals = []

    iters = 0
    convergence = 0
    while solver.successful() and solver.t < sfinal:

        ys = solver.integrate(sfinal, step=True)
        Es, fs, Gs = ravel(ys, ha.n_sp_states)
185
186
        s_vals.append(solver.t)
        E_vals.append(Es)
187
188

        iters += 1
189
190
        # if iters == 176:
        #     break
191
        if iters %10 == 0 and verbose: print("iter: {:>6d} \t scale param: {:0.4f} \t E = {:0.8f}".format(iters, solver.t, Es))
192

193
        if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) < 10**-8 and E_vals[-1] != E_vals[0]:
194
195

            if verbose: print("---- Energy converged at iter {:>06d} with energy {:1.8f}\n".format(iters,E_vals[-1]))
196
197
            convergence = 1
            break
198

199
        if len(E_vals) > 100 and abs(E_vals[-1] - E_vals[-2]) > 1:
200
            if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
201
            break
202

203
        if iters > 20000:
204
            if verbose: print("---- Energy diverged at iter {:>06d} with energy {:3.8f}\n".format(iters,E_vals[-1]))
205
            break
206

207
        if iters % 1000 == 0 and verbose:
208
            print('Iteration {:>06d}'.format(iters))
Jacob Davison's avatar
Jacob Davison committed
209
    flowf = time.time()
210
    end = time.time()
211
212
213
    time_str = "{:2.5f}".format(end-start)

    if verbose: print("IM-SRG(2) converged in {:2.5f} seconds".format(flowf-flowi))
214
215

    del ha, ot, wg, fl, solver, y0, sfinal, ds
Jacob August Davison's avatar
Jacob August Davison committed
216
    
217
    return (convergence, iters, d, g, pb, n_holes+n_particles, s_vals, E_vals, time_str)
218
 
219
220

if __name__ == '__main__':
Jacob Davison's avatar
Jacob Davison committed
221

222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237

    #test_exact('plots_exact_2b/', main)

    main(4,4)

    #holes = int(sys.argv[1])

    # print('convergence,iters,d,g,pb,num states,GS energy,total time')
    # for num_pairs in range(1,21):
    #     #print(n)55
    #     #n = int(num_sp/2)
    #     #holes = num_pairs*2
    #     particles = num_pairs*2

    #     convergence,iters,d,g,pb,sp_states,s_vals,E_vals,time_str = main(holes, particles, verbose=0)
    #     print('{},{},{},{},{},{},{},{}'.format(convergence,iters,d,g,pb,sp_states,E_vals[-1],time_str))
238
    
239
    #     del convergence, iters, d, g, pb, sp_states, s_vals, E_vals, time_str