From 21d49b402a5692ecaf9ff3f0ec736aea7c45f8b6 Mon Sep 17 00:00:00 2001 From: liwt31 Date: Sat, 11 May 2024 11:27:57 +0800 Subject: [PATCH] Update Ohmic SDF in SBM support arbitrary s --- renormalizer/sbm/lib.py | 174 +++++++++++++---------------- renormalizer/sbm/tests/test_sbm.py | 1 + 2 files changed, 77 insertions(+), 98 deletions(-) diff --git a/renormalizer/sbm/lib.py b/renormalizer/sbm/lib.py index f76c8f83..2ed09975 100644 --- a/renormalizer/sbm/lib.py +++ b/renormalizer/sbm/lib.py @@ -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__) @@ -29,7 +29,7 @@ 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 @@ -37,37 +37,65 @@ def func(self, omega_value): 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: @@ -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 @@ -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): """ @@ -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) @@ -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) diff --git a/renormalizer/sbm/tests/test_sbm.py b/renormalizer/sbm/tests/test_sbm.py index 4735f466..c865e3e4 100644 --- a/renormalizer/sbm/tests/test_sbm.py +++ b/renormalizer/sbm/tests/test_sbm.py @@ -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)