Skip to content

Commit

Permalink
Update Ohmic SDF in SBM
Browse files Browse the repository at this point in the history 
support arbitrary s
  • Loading branch information
@liwt31
liwt31 committed May 11, 2024
1 parent 37145a7 commit 21d49b4
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 98 deletions.
174 changes: 76 additions & 98 deletions renormalizer/sbm/lib.py
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
# -*- coding: utf-8 -*-

import logging
import numpy as np
import numpy.polynomial.laguerre as la
import numpy.polynomial.legendre as le
from typing import Union, Tuple

import scipy
import scipy.special
import scipy.optimize

from renormalizer.model import Phonon, SpinBosonModel
from renormalizer.utils import Quantity
from renormalizer.mps.backend import np


logger = logging.getLogger(__name__)
Expand All @@ -29,45 +29,73 @@ def func(self, omega_value):
"""
the function of the Debye-type spectral density function
"""
return 2. * self.lamb * omega_value * self.omega_c / (omega_value ** 2 + self.omega_c ** 2)
return 2.0 * self.lamb * omega_value * self.omega_c / (omega_value**2 + self.omega_c**2)


DebyeSDF = DebyeSpectralDensityFunction


class SpectralDensityFunction:
r"""
the ohmic spectral density function
The Ohmic spectral density function
J(\omega) = \pi / 2 \alpha \omega e^{-\omega/\omega_c}
"""

def __init__(self, alpha, omega_c: Quantity):
def __init__(self, alpha: float, omega_c: Union[Quantity, float], s: float = 1):
self.alpha = alpha
self.omega_c = omega_c.as_au()
if isinstance(omega_c, Quantity):
omega_c = omega_c.as_au()
self.omega_c = omega_c
self.s = s

def reno(self, omega_l) -> float:
def integrate_func(x):
return self.func(x) / x**2

res = scipy.integrate.quad(integrate_func, a=omega_l, b=self.omega_c * 30)
logger.info(f"integrate: {res[0]}, {res[1]}")
re = np.exp(-res[0] * 2 / np.pi)

return re

def adiabatic_renormalization(self, delta_quan: Quantity, p: float):
delta = delta_quan.as_au()
def adiabatic_renormalization(self, delta: Union[Quantity, float], p: float) -> Tuple[float, float]:
"""
the cut-off omega_l is p*delta
"""
if isinstance(delta, Quantity):
delta = delta.as_au()
loop = 0
re = 1.
re = 1.0
while loop < 50:
re_old = re
omega_l = delta * re * p
re = np.exp(-self.alpha * scipy.special.expn(1, omega_l / self.omega_c))

def integrate_func(x):
return self.func(x) / x**2

res = scipy.integrate.quad(integrate_func, a=omega_l, b=self.omega_c * 30)
logger.info(f"integrate: {res[0]}, {res[1]}")
re = np.exp(-res[0] * 2 / np.pi)
loop += 1
logger.info(f"re, {re_old}, {re}")
if np.allclose(re, re_old):
break

return delta_quan * re, Quantity(delta * re * p)
return delta * re, delta * re * p

def func(self, omega_value):
"""
the function of the ohmic spectral density function
"""
return np.pi / 2. * self.alpha * omega_value * np.exp(-omega_value / self.omega_c)
return (
np.pi / 2.0 * self.alpha
* omega_value**self.s * self.omega_c ** (1 - self.s)
* np.exp(-omega_value / self.omega_c)
)

@staticmethod
def post_process(omega_value, c_j2, ifsort):
displacement_array = np.sqrt(c_j2) / omega_value ** 2
def post_process(omega_value, c_j2, ifsort=True):
displacement_array = np.sqrt(c_j2) / omega_value**2
if ifsort:
idx = np.argsort(c_j2 / omega_value)[::-1]
else:
Expand All @@ -85,80 +113,27 @@ def _dos_Wang1(self, nb, omega_value):
"""
return (nb + 1) / self.omega_c * np.exp(-omega_value / self.omega_c)

def Wang1(self, nb, ifsort=True):
"""
def Wang1(self, nb):
r"""
Wang's 1st scheme discretization
"""
omega_value = np.array([-np.log(-float(j) / (nb + 1) + 1.) * self.omega_c for j in range(1, nb + 1, 1)])
omega_value = np.array([-np.log(-float(j) / (nb + 1) + 1.0) * self.omega_c for j in range(1, nb + 1, 1)])

# general form
# c_j2 = 2./np.pi * omega_value * self.func(omega_value) / self._dos_Wang1(nb, omega_value)

# excat form
c_j2 = omega_value ** 2 * self.alpha * self.omega_c / (nb + 1)

return self.post_process(omega_value, c_j2, ifsort)

def legendre(self, nb, x0, x1, ifsort=True):
"""
Legendre polynomial fit [x0, x1] to [-1,1]
omega_m is the cutoff
"""
omega_value, w = le.leggauss(nb)
omega_value = (omega_value + (x1 + x0) / (x1 - x0)) * (x1 - x0) / 2.
c_j2 = w * (x1 - x0) / 2. * self.alpha * omega_value ** 2 * np.exp(-omega_value / self.omega_c)

return self.post_process(omega_value, c_j2, ifsort)

def laguerre(self, nb, ifsort=True):
assert nb <= 100
c_j2 = 2.0 / np.pi * omega_value * self.func(omega_value) / self._dos_Wang1(nb, omega_value)

omega_value, w = la.laggauss(nb)
omega_value *= self.omega_c
c_j2 = w * self.alpha * self.omega_c * omega_value ** 2

return self.post_process(omega_value, c_j2, ifsort)
return omega_value, c_j2

def trapz(self, nb, x0, x1, ifsort=True):
def trapz(self, nb, x0, x1):
dw = (x1 - x0) / float(nb)
xlist = [x0 + i * dw for i in range(nb + 1)]
omega_value = np.array([(xlist[i] + xlist[i + 1]) / 2. for i in range(nb)])
c_j2 = np.array([(self.func(xlist[i]) + self.func(xlist[i + 1])) / 2 for i in range(nb)]) * 2. / np.pi * omega_value * dw

return self.post_process(omega_value, c_j2, ifsort)

def _opt_cut(self, p):
"""
p is the percent of the height of the SDF
"""
assert 0 < p < 1
cut = np.exp(-1.) * p
omega_value = np.array([(xlist[i] + xlist[i + 1]) / 2.0 for i in range(nb)])
c_j2 = (
np.array([(self.func(xlist[i]) + self.func(xlist[i + 1])) / 2 for i in range(nb)])
* 2.0 / np.pi * omega_value * dw
)

def F(x):
return x * np.exp(-x) - cut

def fprime(x):
return np.exp(-x) * (1. - x)

def fprime2(x):
return (x - 2.) * np.exp(-x)

x1 = scipy.optimize.newton(F, 0.0, fprime=fprime, fprime2=fprime2)
x2 = scipy.optimize.newton(F, -np.log(cut), fprime=fprime, fprime2=fprime2)

return x1 * self.omega_c, x2 * self.omega_c

def plot_data(self, x0, x1, n, omega_value, c_j2, sigma=0.1):
"""
plot the spectral density function (continuous and discrete)
"""
x = np.linspace(x0, x1, n)
y_c = self.func(x)

y_d = np.einsum("i,ji -> j", c_j2 / omega_value * np.pi / 2. * 1 / np.sqrt(2 * np.pi * sigma ** 2),
np.exp(-(np.subtract.outer(x, omega_value) / sigma) ** 2 / 2.))

return x, y_c, y_d
return omega_value, c_j2


OhmicSDF = SpectralDensityFunction
Expand All @@ -175,26 +150,22 @@ def __init__(self, ita, omega_c, beta, omega_limit):
self.beta = beta
self.omega_limit = omega_limit


def reno(self, omega_l):
def integrate_func(x):
return self.func(x) / x**2

res = scipy.integrate.quad(integrate_func, a=omega_l,
b=omega_l*1000)
res = scipy.integrate.quad(integrate_func, a=omega_l, b=omega_l * 1000)
logger.info(f"integrate: {res[0]}, {res[1]}")
re = np.exp(-res[0]*2/np.pi)
re = np.exp(-res[0] * 2 / np.pi)

return re


def func(self, omega_value):
"""
the function of the spectral density function
"""
theta = np.arctan(omega_value/self.omega_c)
return self.ita * np.sin(self.beta * theta) / (1 + omega_value**2/self.omega_c**2) ** (self.beta / 2)

theta = np.arctan(omega_value / self.omega_c)
return self.ita * np.sin(self.beta * theta) / (1 + omega_value**2 / self.omega_c**2) ** (self.beta / 2)

def _dos_Wang1(self, A, omega_value):
"""
Expand All @@ -207,11 +178,13 @@ def Wang1(self, nb):
"""
Wang's 1st scheme discretization
"""

def integrate_func(x):
return self.func(x) / x
A = (nb + 1 ) / scipy.integrate.quad(integrate_func, a=0, b=self.omega_limit)[0]

A = (nb + 1) / scipy.integrate.quad(integrate_func, a=0, b=self.omega_limit)[0]
logger.info(scipy.integrate.quad(integrate_func, a=0, b=self.omega_limit)[0] * 4 / np.pi)
logger.info(2*self.ita)
logger.info(2 * self.ita)
nsamples = int(1e7)
delta = self.omega_limit / nsamples
omega_value_big = np.linspace(delta, self.omega_limit, nsamples)
Expand All @@ -224,16 +197,21 @@ def integrate_func(x):
assert len(omega_value) == nb

# general form
c_j2 = 2./np.pi * omega_value * self.func(omega_value) / self._dos_Wang1(A, omega_value)

c_j2 = 2.0 / np.pi * omega_value * self.func(omega_value) / self._dos_Wang1(A, omega_value)

return omega_value, c_j2


def param2mollist(alpha: float, raw_delta: Quantity, omega_c: Quantity, renormalization_p: float, n_phonons: int):
sdf = SpectralDensityFunction(alpha, omega_c)
def param2mollist(
alpha: float,
raw_delta: Quantity,
omega_c: Quantity,
renormalization_p: float,
n_phonons: int,
):
sdf = SpectralDensityFunction(alpha, omega_c, s=1)
delta, max_omega = sdf.adiabatic_renormalization(raw_delta, renormalization_p)
omega_list, displacement_list = sdf.trapz(n_phonons, 0.0, max_omega.as_au())

ph_list = [Phonon.simplest_phonon(o, d) for o,d in zip(omega_list, displacement_list)]
return SpinBosonModel(Quantity(0), delta, ph_list)
omega_list, displacement_list = sdf.trapz(n_phonons, 0.0, max_omega)
omega_list, displacement_list = sdf.post_process(omega_list, displacement_list)
ph_list = [Phonon.simplest_phonon(o, d) for o, d in zip(omega_list, displacement_list)]
return SpinBosonModel(Quantity(0), Quantity(delta), ph_list)
1 change: 1 addition & 0 deletions renormalizer/sbm/tests/test_sbm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ def test_sdf():
omega_c = Quantity(5)
sdf = SpectralDensityFunction(alpha, omega_c)
omega_list, displacement_list = sdf.trapz(200, 0.0, 50)
omega_list, displacement_list = sdf.post_process(omega_list, displacement_list)

ph_list = [Phonon.simplest_phonon(o, d) for o,d in zip(omega_list, displacement_list)]
mol_reor = sum(ph.reorganization_energy.as_au() for ph in ph_list)
Expand Down

0 comments on commit 21d49b4

Please sign in to comment.