lime/models/pyrazine.py
def polariton_hamiltonian(n_el, n_vc, n_vt, cav, g):
    """
    contruct vibronic-cavity basis for polaritonic states, that is a
    direct product of electron-vibration-photon space
    """

    #n_basis = n_el * n_cav * n_vc * n_vt

    freq_vc = 952. * wavenumber2hartree
    freq_vt = 597. * wavenumber2hartree

    Eshift = np.array([0.0, 31800.0, 39000]) * wavenumber2hartree
    kappa = np.array([0.0, -847.0, 1202.]) * wavenumber2hartree
    coup = 2110.0 * wavenumber2hartree # inter-state coupling lambda


    # indentity matrices in each subspace
    I_el = identity(n_el)
    I_cav =  cav.idm
    I_vc = identity(n_vc)
    I_vt = identity(n_vt)

    h_cav = cav.ham()
    h_vt = boson(freq_vt, n_vt, ZPE=False)
    h_vc = boson(freq_vc, n_vc, ZPE=False)

    h_el = np.diagflat(Eshift)

    # the bare term in the system Hamiltonian
    h0 = kron(h_el, kron(I_cav, kron(I_vc, I_vt))) + \
        kron(I_el, kron(h_cav, kron(I_vc, I_vt))) \
        + kron(I_el, kron(I_cav, kron(h_vc, I_vt))) +\
        kron(I_el, kron(I_cav, kron(I_vc, h_vt))) \

    X = pos(n_vt)


    h1 = kron(np.diagflat(kappa), kron(I_cav, kron(I_vc, X)))


    Xc = pos(n_vc)

    trans_el = np.zeros((n_el, n_el)) # electronic excitation operator
    #deex = np.zeros((n_el, n_el)) # deexcitation
    #deex[2, 1] = 1.0
    #ex[1, 2] = 1.0
    trans_el[1,2] = trans_el[2,1] = 1.0 #= ex + deex

    #h_fake = kron(np.diagflat(kappa), kron(I_cav, kron(Xc, I_vt)))

    ###
#    h_m = np.zeros((n_basis, n_basis))
#
#    for m, b0 in enumerate(basis_set):
#        for n, b1 in enumerate(basis_set):
##            h_m[m,n] = h_el[b0.n_el, b1.n_el] * h_cav[b0.n_cav, b1.n_cav] * h_vc[b0.n_vc, b1.n_vc]
#            h_m[m,n] = trans_el[b0.n_el, b1.n_el] * I_cav[b0.n_cav, b1.n_cav] \
#                    * Xc[b0.n_vc, b1.n_vc] * I_vt[b0.n_vt, b1.n_vt]

    h2 = coup * kron(trans_el, kron(I_cav, kron(Xc, I_vt)), format='csr')

    # if n_cav = n _el
    deex_cav = cav.get_annihilate()
    ex_cav = cav.get_create()

    d_ex = np.zeros((n_el, n_el)) # electronic excitation operator
    d_deex = np.zeros((n_el, n_el)) # deexcitation
    d_deex[0, 2] = 1.0
    d_ex[2, 0] = 1.0

    dip = d_deex + d_ex

    h3 = g * kron(dip, kron(deex_cav + ex_cav, kron(I_vc, I_vt)))

    h_s = h0 + h1 + h2 + h3

    # polaritonic states can be obtained by diagonalizing H_S
    # v is the basis transformation matrix, v[i,j] = <old basis i| polaritonic state j>
    #eigvals, v = np.linalg.eigh(h_s)

    #h_s = csr_matrix(h_s)

    # collapse operators in dissipative dynamics
    Sc = kron(I_el, kron(I_cav, kron(Xc, I_vt)), format='csr')
    St = kron(I_el, kron(I_cav, kron(I_vc, X)), format='csr')

#    St = csr_matrix(St)
#    Sc = csr_matrix(Sc)

    return h_s, Sc, St

    contruct vibronic-cavity basis for polaritonic states, that is a
    direct product of electron-vibration-photon space
    
contruct vibronic-cavity basis for polaritonic states, that is a
direct product of electron-vibration-photon space
contruct vibronic - cavity basis for polaritonic states that is a direct product of electron - vibration - photon space
def polariton_hamiltonian(n_el, n_vc, n_vt, cav, g):
    freq_vc = 952. * wavenumber2hartree
    freq_vt = 597. * wavenumber2hartree
    Eshift = np.array([0.0, 31800.0, 39000]) * wavenumber2hartree
    kappa = np.array([0.0, -847.0, 1202.]) * wavenumber2hartree
    coup = 2110.0 * wavenumber2hartree 
    I_el = identity(n_el)
    I_cav =  cav.idm
    I_vc = identity(n_vc)
    I_vt = identity(n_vt)
    h_cav = cav.ham()
    h_vt = boson(freq_vt, n_vt, ZPE=False)
    h_vc = boson(freq_vc, n_vc, ZPE=False)
    h_el = np.diagflat(Eshift)
    h0 = kron(h_el, kron(I_cav, kron(I_vc, I_vt))) + \
        kron(I_el, kron(h_cav, kron(I_vc, I_vt))) \
        + kron(I_el, kron(I_cav, kron(h_vc, I_vt))) +\
        kron(I_el, kron(I_cav, kron(I_vc, h_vt))) \
    X = pos(n_vt)
    h1 = kron(np.diagflat(kappa), kron(I_cav, kron(I_vc, X)))
    Xc = pos(n_vc)
    trans_el = np.zeros((n_el, n_el)) 
    trans_el[1,2] = trans_el[2,1] = 1.0 
    for m, b0 in enumerate(basis_set):
        for n, b1 in enumerate(basis_set):
            h_m[m,n] = h_el[b0.n_el, b1.n_el] * h_cav[b0.n_cav, b1.n_cav] * h_vc[b0.n_vc, b1.n_vc]
            h_m[m,n] = trans_el[b0.n_el, b1.n_el] * I_cav[b0.n_cav, b1.n_cav] \
                    * Xc[b0.n_vc, b1.n_vc] * I_vt[b0.n_vt, b1.n_vt]
    h2 = coup * kron(trans_el, kron(I_cav, kron(Xc, I_vt)), format='csr')
     if n_cav = n _el
    deex_cav = cav.get_annihilate()
    ex_cav = cav.get_create()
    d_ex = np.zeros((n_el, n_el)) 
    d_deex = np.zeros((n_el, n_el))  deexcitation
    d_deex[0, 2] = 1.0
    d_ex[2, 0] = 1.0
    dip = d_deex + d_ex
    h3 = g * kron(dip, kron(deex_cav + ex_cav, kron(I_vc, I_vt)))
    h_s = h0 + h1 + h2 + h3
     polaritonic states can be obtained by diagonalizing H_S
     v is the basis transformation matrix, v[i,j] = <old basis i| polaritonic state j>
    eigvals, v = np.linalg.eigh(h_s)
    h_s = csr_matrix(h_s)
     collapse operators in dissipative dynamics
    Sc = kron(I_el, kron(I_cav, kron(Xc, I_vt)), format='csr')
    St = kron(I_el, kron(I_cav, kron(I_vc, X)), format='csr')
    St = csr_matrix(St)
    Sc = csr_matrix(Sc)
    return h_s, Sc, St
def polariton_hamiltonian ( n_el , n_vc , n_vt , cav , g ) : freq_vc = 952. * wavenumber2hartree freq_vt = 597. * wavenumber2hartree Eshift = np . array ( [ 0.0 , 31800.0 , 39000 ] ) * wavenumber2hartree kappa = np . array ( [ 0.0 , - 847.0 , 1202. ] ) * wavenumber2hartree coup = 2110.0 * wavenumber2hartree I_el = identity ( n_el ) I_cav = cav . idm I_vc = identity ( n_vc ) I_vt = identity ( n_vt ) h_cav = cav . ham ( ) h_vt = boson ( freq_vt , n_vt , ZPE = False ) h_vc = boson ( freq_vc , n_vc , ZPE = False ) h_el = np . diagflat ( Eshift ) h0 = kron ( h_el , kron ( I_cav , kron ( I_vc , I_vt ) ) ) + kron ( I_el , kron ( h_cav , kron ( I_vc , I_vt ) ) ) + kron ( I_el , kron ( I_cav , kron ( h_vc , I_vt ) ) ) + kron ( I_el , kron ( I_cav , kron ( I_vc , h_vt ) ) ) X = pos ( n_vt ) h1 = kron ( np . diagflat ( kappa ) , kron ( I_cav , kron ( I_vc , X ) ) ) Xc = pos ( n_vc ) trans_el = np . zeros ( ( n_el , n_el ) ) trans_el [ 1 , 2 ] = trans_el [ 2 , 1 ] = 1.0 h2 = coup * kron ( trans_el , kron ( I_cav , kron ( Xc , I_vt ) ) , format = 'csr' ) deex_cav = cav . get_annihilate ( ) ex_cav = cav . get_create ( ) d_ex = np . zeros ( ( n_el , n_el ) ) d_deex = np . zeros ( ( n_el , n_el ) ) d_deex [ 0 , 2 ] = 1.0 d_ex [ 2 , 0 ] = 1.0 dip = d_deex + d_ex h3 = g * kron ( dip , kron ( deex_cav + ex_cav , kron ( I_vc , I_vt ) ) ) h_s = h0 + h1 + h2 + h3 Sc = kron ( I_el , kron ( I_cav , kron ( Xc , I_vt ) ) , format = 'csr' ) St = kron ( I_el , kron ( I_cav , kron ( I_vc , X ) ) , format = 'csr' ) return h_s , Sc , St
contruct vibronic-cavity basis for polaritonic states, that is a
direct product of electron-vibration-photon space
contruct vibronic - cavity basis for polaritonic states that is a direct product of electron - vibration - photon space
"""
    contruct vibronic-cavity basis for polaritonic states, that is a
    direct product of electron-vibration-photon space
    """ #n_basis = n_el * n_cav * n_vc * n_vt # inter-state coupling lambda # indentity matrices in each subspace # the bare term in the system Hamiltonian # electronic excitation operator #deex = np.zeros((n_el, n_el)) # deexcitation #deex[2, 1] = 1.0 #ex[1, 2] = 1.0 #= ex + deex #h_fake = kron(np.diagflat(kappa), kron(I_cav, kron(Xc, I_vt))) ### #    h_m = np.zeros((n_basis, n_basis)) # #    for m, b0 in enumerate(basis_set): #        for n, b1 in enumerate(basis_set): ##            h_m[m,n] = h_el[b0.n_el, b1.n_el] * h_cav[b0.n_cav, b1.n_cav] * h_vc[b0.n_vc, b1.n_vc] #            h_m[m,n] = trans_el[b0.n_el, b1.n_el] * I_cav[b0.n_cav, b1.n_cav] \ #                    * Xc[b0.n_vc, b1.n_vc] * I_vt[b0.n_vt, b1.n_vt] # if n_cav = n _el # electronic excitation operator # deexcitation # polaritonic states can be obtained by diagonalizing H_S # v is the basis transformation matrix, v[i,j] = <old basis i| polaritonic state j> #eigvals, v = np.linalg.eigh(h_s) #h_s = csr_matrix(h_s) # collapse operators in dissipative dynamics #    St = csr_matrix(St) #    Sc = csr_matrix(Sc)
n_el
n_vc
n_vt
cav
g