added generalised functions to compute dynamics of multiple systems (…

…with a field)

fixed contraction.py, generalised dynamics.py for super-system, added examples

removed unnecessary csv

scripts for two-time correlations - needs testing!

remove old testing script

Fix setup.py

Version bump to 0.3.1

added disordered mean field code

added code to check photon number convergence for large system numbers

examples/mean-field_correlations/dynamics.py

@@ -0,0 +1,238 @@
+#!/usr/bin/env python
+
+"""Rough script to measure two-time correlators for mean-field Hamiltonian"""
+
+import pickle, os, sys
+sys.path.insert(0,'../../') # Make OQuPy accessible
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import oqupy
+import oqupy.operators as op
+
+
+# Output files
+data_dir = 'data'
+fig_dir = 'figures'
+dynamics_plotfp = os.path.join(fig_dir, 'dynamics_plot.pdf')
+datafp = os.path.join(data_dir, 'dynamics_correlators.pkl')
+if not os.path.isdir(data_dir):
+    os.makedirs(data_dir)
+if not os.path.isdir(fig_dir):
+    os.makedirs(fig_dir)
+
+
+# Spin operators
+sigma_z = oqupy.operators.sigma('z')
+sigma_p = oqupy.operators.sigma('+')
+sigma_m = oqupy.operators.sigma('-')
+
+# Bath parameters
+nu_c = 0.15
+a = 0.25
+T = 0.026
+
+# System parameters
+dim = 2
+wc = 0.0
+gn = 0.2    
+gam_down = 0.01 
+gam_up   = 0.01
+kappa = 0.01
+
+w0 = 0.0
+
+# Initial conditions
+initial_state = oqupy.operators.spin_dm('z-')
+initial_field = np.sqrt(0.05)
+
+# Computational parameters
+# N.B. Not sensible values!
+ts = 400 # steady-state time
+tp = 100 # time to measure field period from (should really be in steady-state by tp)
+tf = 600 # final time
+rotating_frame_freq = None # record rotating frame freq. (used below)
+dt = 0.5
+dkmax = 20
+epsrel = 10**(-4)
+
+
+local_times = [0.0]
+local_fields = [initial_field]
+def store_local_field(t, a):
+    if not np.isclose(t, local_times[-1]+dt):
+        # in input parsing random times are passed to field_eom, avoid recording
+        # these
+        return
+    if t > local_times[-1]:
+        local_times.append(t)
+        local_fields.append(a)
+
+# Measure period of oscillation and adjust global variables w0 and wc
+# For multiple systems would need to adjust each w0
+def move_to_rotating_frame():
+    global rotating_frame_freq, w0, wc
+    start_step = round(tp/dt)
+    end_step = None
+    field_samples = local_fields[start_step:end_step]
+    # Array of periods measured in steps, taking each period as the time between 3 intercepts of the horizontal axis
+    period_steps = []
+    # count number of intercepts
+    intercepts = 0
+    # on first intercept, record number of steps
+    recorded_step = 0
+    # determine where sign of real part of field changes (assume evolution continuous)
+    sign_changes = np.diff(np.sign(np.real(field_samples)))
+    for step, change in enumerate(sign_changes):
+        # If sign changes, we have an intercept
+        if change != 0:
+            intercepts+=1
+            # record step of first intercept (3 intercepts make 1 period)
+            if intercepts==1:
+                recorded_step=step
+        if intercepts==3:
+        # Period is difference between step of third intercept and step of first intercept
+            period_steps.append(step-recorded_step)
+            # reset counter; hopefully measure multiple periods and average to minimise numerical error
+            # due to timestep not exactly aligning with intercepts
+            intercepts=0
+    num_periods = len(period_steps)
+    if num_periods == 0:
+        # Nothing to do; no periods measured (field not oscillatory)
+        print('\nNo field oscillations recorded between t={} and t={}'.format(tp, ts)) 
+        rotating_frame_freq = 0.0
+        return
+    elif num_periods <= 5:
+        print('\nOnly {} periods recorded between t={} and t={} - rotating frame '\
+                'frequency may be inaccurate.'.format(num_periods, tp, ts))
+    # average period in units time (not steps)
+    average_period = dt * np.average(period_steps)
+    lasing_angular_freq = 2*np.pi / average_period
+    phi0 = np.angle(field_samples[-2])
+    phi1 = np.angle(field_samples[-1])
+    # whether phase is increasing or decreasing
+    lasing_direction = np.sign(phi1-phi0)
+    # np.angle has discontinuity on negative Im axis, so above fails if phi0 in upper left quadrant and phi1 in bottom left
+    if phi1<-np.pi/2 and phi0>np.pi/2:
+        lasing_direction = -1
+    # add corresponding angular frequency from both w0 and wc. This should result in a stationary solution
+    # (add as alpha rotates at negative of rotating frame freq)
+    rotating_frame_freq = lasing_direction*lasing_angular_freq
+    w0 += rotating_frame_freq  # MUTLI-SYSTEM GENERALISATION ?
+    wc += rotating_frame_freq
+    print('Adjusted w0 and wc by rotating_frame_freq {:.3f}'.format(rotating_frame_freq))
+
+
+# Functions passed to tempo
+def field_eom(t, state, a):
+    if rotating_frame_freq is None:
+        # need to keep a local copy of field so can calculate period of
+        # oscillation later. TODO: implement in OQuPy itself 
+        store_local_field(t, a)
+    if t >= ts:
+        # In steady-state, move to rotating frame and stop evolving field
+        if rotating_frame_freq is None:
+            move_to_rotating_frame()
+        return 0.0
+    expect_val = np.matmul(sigma_m, state).trace()
+    return -(1j * wc + kappa) * a - 0.5j * gn * expect_val
+def H_MF(t, a):
+    return 0.5 * w0 * sigma_z +\
+        0.5 * gn * (a * sigma_p + np.conj(a) * sigma_m)
+
+
+system = oqupy.TimeDependentSystemWithField(H_MF,
+        field_eom,
+        gammas = [lambda t: gam_down, lambda t: gam_up],
+        lindblad_operators = [lambda t: sigma_m, lambda t: sigma_p]
+        )
+correlations = oqupy.PowerLawSD(alpha=a,
+                                zeta=1,
+                                cutoff=nu_c,
+                                cutoff_type='gaussian',
+                                temperature=T)
+bath = oqupy.Bath(0.5 * sigma_z, correlations)
+
+pt_tempo_parameters = oqupy.TempoParameters(
+    dt=dt,
+    dkmax=dkmax,
+    epsrel=epsrel)
+
+# compute PT to final time tf
+process_tensor = oqupy.pt_tempo_compute(bath=bath,
+                                        start_time=0.0,
+                                        end_time=tf,
+                                        parameters=pt_tempo_parameters)
+
+# Control objects for two-time correlator measurement
+control_sm = oqupy.Control(dim)
+control_sp = oqupy.Control(dim)
+# N.B. ts must be a float otherwise (if int) interpreted as timestep 
+control_sm.add_single(float(ts), op.left_super(sigma_m))
+control_sp.add_single(float(ts), op.left_super(sigma_p))
+
+# Two sets of dynamics, one for each two-time correlator
+dynamics_sm = oqupy.compute_dynamics_with_field(
+    system=system,
+    process_tensor=process_tensor,
+	initial_field=initial_field,
+    control=control_sm,
+    start_time=0.0,
+    initial_state=initial_state)
+times, sp = dynamics_sm.expectations(oqupy.operators.sigma('+')/2, real=False)
+ts_index = next((i for i, t in enumerate(times) if t >= ts), None)
+corr_times = times[ts_index:] - ts
+spsm = sp[ts_index:] 
+first_rotating_frame_freq = rotating_frame_freq
+# reset rotating frame frequency and local storage variables
+rotating_frame_freq = None 
+local_times = [0.0]
+local_fields = [initial_field]
+dynamics_sp = oqupy.compute_dynamics_with_field(
+    system=system,
+    process_tensor=process_tensor,
+	initial_field=initial_field,
+    control=control_sp,
+    start_time=0.0,
+    initial_state=initial_state)
+times, sm = dynamics_sp.expectations(oqupy.operators.sigma('-')/2, real=False)
+smsp = sm[ts_index:] # <sigma^-(t) sigma^+(0)>
+# consistency check
+assert rotating_frame_freq == first_rotating_frame_freq
+assert len(smsp) == len(spsm) == len(corr_times)
+
+# save truncated times, correlators and parameters used by spectrum.py
+save_dic = {
+        'times': corr_times,
+        'spsm': spsm,
+        'smsp': smsp,
+        'params': {
+            'dt': dt,
+            'wc': wc,
+            'w0': w0,
+            'rotating_frame_freq': rotating_frame_freq,
+            'kappa': kappa,
+            'gn': gn,
+            }
+        }
+# 
+with open(datafp, 'wb') as fb:
+    pickle.dump(save_dic, fb)
+print('Times and correlator values saved to {}'.format(datafp))
+
+# Plot fields and polarisation
+times, s_z = dynamics_sm.expectations(oqupy.operators.sigma('z')/2, real=True)
+times, fields = dynamics_sm.field_expectations()
+fig, axes = plt.subplots(2, figsize=(9,10))
+axes[0].plot(times, s_z)
+axes[1].plot(times, np.real(fields))
+axes[1].set_xlabel('t')
+axes[0].set_ylabel('<Sz>')
+axes[1].set_ylabel('Re<a>')
+axes[1].axvline(x=tp, c='g') # corresponds to time measure period from
+axes[1].axvline(x=ts, c='r') # corresponds to time measure correlators from
+fig.savefig(dynamics_plotfp, bbox_inches='tight')
+
+
+

examples/mean-field_correlations/spectrum.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python
+
+"""Script to calculate spectral weight and photoluminescence from two-time correlators"""
+
+import pickle, sys
+
+import numpy as np
+import matplotlib.pyplot as plt
+plt.rc('text', usetex=True)
+plt.rc('font', **{'size':14})
+plt.rc('text.latex', preamble=r'\usepackage{amsmath}')
+
+# input file
+inputfp = 'data/dynamics_correlators.pkl'
+
+# load
+with open(inputfp, 'rb') as fb:
+    data = pickle.load(fb)
+times = data['times']
+spsm = data['spsm']
+smsp = data['smsp']
+params = data['params']
+
+dt = params['dt']
+# calculations in rotating frame, plot in original frame
+wc_rotating = params['wc'] + params['rotating_frame_freq']
+w0_rotating = params['w0'] + params['rotating_frame_freq']
+nus = 2 * np.pi * np.fft.fftshift(np.fft.fftfreq(len(times), d=dt))
+plot_nus = 1e3 * (nus - params['rotating_frame_freq']) # units meV
+
+
+# calculate self-energies
+fft_smsp = np.fft.fftshift(dt * np.fft.ifft(smsp, norm='forward'))
+fft_spsm_conjugate = np.fft.fftshift(dt * np.fft.ifft(np.conjugate(spsm), norm='forward'))
+# Sigma^{-+}
+energy_mp = -(1j/4)*params['gn']**2*(fft_smsp-fft_spsm_conjugate)
+# Sigma^{--} 
+energy_mm = -(1j/2)*params['gn']**2*(np.real(fft_smsp+fft_spsm_conjugate))
+
+
+# calculate unperturbed Green's functions
+# inverse of non-interacting retarded function
+def D0RI(nu):
+    return nu-wc_rotating+1j*params['kappa']
+# N.B. DOIK = - DORI * DOK * DOA happens to be a constant
+def D0IK(nu):
+    return 2j * params['kappa']
+
+# calculate interacting green's functions
+inverse_retarded = D0RI(nus) - energy_mp
+keldysh_inverse = D0IK(nus) - energy_mm
+retarded = 1/inverse_retarded
+advanced = np.conjugate(retarded)
+keldysh = - retarded * keldysh_inverse * advanced
+
+# calculate pl and spectral weight
+#pl = (1j/2) * (keldysh - retarded + advanced)
+# use explicit formula for PL in terms of self energies 
+pl = (np.imag(energy_mp) - 0.5*np.imag(energy_mm)) \
+        / np.abs((nus - wc_rotating) + 1j * params['kappa'] - energy_mp)**2
+# spectral weight
+spectral_weight = -2*np.imag(retarded)
+
+
+# Plot correlators, spectral weight, photoluminescence 
+fig1, axes1 = plt.subplots(2, figsize=(9,6), sharex=True)
+fig2, ax2 = plt.subplots(figsize=(9,4))
+fig3, ax3 = plt.subplots(figsize=(9,4))
+
+axes1[1].set_xlabel('\(t\)')
+axes1[0].plot(times, np.real(smsp), label=r'\(\text{Re}\langle \sigma^-(t) \sigma^+(0) \rangle\)')
+axes1[0].plot(times, np.imag(smsp), label=r'\(\text{Im}\langle \sigma^-(t) \sigma^+(0) \rangle\)')
+axes1[1].plot(times, np.real(spsm), label=r'\(\text{Re}\langle \sigma^+(t) \sigma^-(0) \rangle\)')
+axes1[1].plot(times, np.imag(spsm), label=r'\(\text{Im}\langle \sigma^+(t) \sigma^-(0) \rangle\)')
+axes1[0].legend()
+axes1[1].legend()
+ax2.set_xlabel(r'\(\nu\)')
+ax2.set_ylabel(r'\(\varrho\)', rotation=0, labelpad=20)
+ax2.plot(plot_nus, spectral_weight)
+ax2.set_xlim([-300,300])
+ax3.set_xlabel(r'\(\nu\)')
+ax3.set_ylabel(r'\(\mathcal{L}\)', rotation=0, labelpad=20)
+ax3.plot(plot_nus, pl)
+ax3.set_xlim([-300,300])
+fig1.savefig('figures/correlators.pdf', bbox_inches='tight')
+fig2.savefig('figures/spectral_weight.pdf', bbox_inches='tight')
+fig3.savefig('figures/photoluminescence.pdf', bbox_inches='tight')