diff --git a/src/external/BMWtools/BMWIntegrator.cpp b/src/external/BMWtools/BMWIntegrator.cpp index 96ce4154..bd810c42 100644 --- a/src/external/BMWtools/BMWIntegrator.cpp +++ b/src/external/BMWtools/BMWIntegrator.cpp @@ -35,6 +35,182 @@ #define SEED 0 #define STATEFILE NULL +//----------------------------------------------------------------------------- + +std::vector TPointPWaveGapIntegralCuhre::fPar; + +//----------------------------------------------------------------------------- +/** + *

Integrate the function using the Cuhre interface + * + *

return: + * - value of the integral + */ +double TPointPWaveGapIntegralCuhre::IntegrateFunc(int tag) +{ + const unsigned int NCOMP(1); + const unsigned int NVEC(1); + const double EPSREL (1e-4); + const double EPSABS (1e-6); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (0); + const unsigned int MAXEVAL (50000); + + const unsigned int KEY (13); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + if (tag == 0) + Cuhre(fNDim, NCOMP, Integrand_aa, USERDATA, NVEC, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, STATEFILE, SPIN, + &nregions, &neval, &fail, integral, error, prob); + else + Cuhre(fNDim, NCOMP, Integrand_cc, USERDATA, NVEC, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, STATEFILE, SPIN, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave point, aa==bb component + * + *

return: + * - 0 + * + * \param ndim number of dimensions of the integral (2 here) + * \param x point where the function should be evaluated + * \param ncomp number of components of the integrand (1 here) + * \param f function value + * \param userdata additional user parameters (required by the interface, NULL here) + */ +int TPointPWaveGapIntegralCuhre::Integrand_aa(const int *ndim, const double x[], + const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc} +{ + double z = x[1]*fPar[3]; + double deltasq(pow(sqrt(1.0-z*z)*fPar[1],2.0)); + f[0] = (1.0-z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + return 0; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave point, cc component + * + *

return: + * - 0 + * + * \param ndim number of dimensions of the integral (2 here) + * \param x point where the function should be evaluated + * \param ncomp number of components of the integrand (1 here) + * \param f function value + * \param userdata additional user parameters (required by the interface, NULL here) + */ +int TPointPWaveGapIntegralCuhre::Integrand_cc(const int *ndim, const double x[], + const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc} +{ + double z = x[1]*fPar[3]; + double deltasq(pow(sqrt(1.0-z*z)*fPar[1],2.0)); + f[0] = (z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + return 0; +} + +//----------------------------------------------------------------------------- + +std::vector TLinePWaveGapIntegralCuhre::fPar; + +//----------------------------------------------------------------------------- +/** + *

Integrate the function using the Cuhre interface + * + *

return: + * - value of the integral + */ +double TLinePWaveGapIntegralCuhre::IntegrateFunc(int tag) +{ + const unsigned int NCOMP(1); + const unsigned int NVEC(1); + const double EPSREL (1e-4); + const double EPSABS (1e-6); + const unsigned int VERBOSE (0); + const unsigned int LAST (4); + const unsigned int MINEVAL (0); + const unsigned int MAXEVAL (50000); + + const unsigned int KEY (13); + + int nregions, neval, fail; + double integral[NCOMP], error[NCOMP], prob[NCOMP]; + + if (tag == 0) + Cuhre(fNDim, NCOMP, Integrand_aa, USERDATA, NVEC, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, STATEFILE, SPIN, + &nregions, &neval, &fail, integral, error, prob); + else + Cuhre(fNDim, NCOMP, Integrand_cc, USERDATA, NVEC, + EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, + KEY, STATEFILE, SPIN, + &nregions, &neval, &fail, integral, error, prob); + + return integral[0]; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave line, aa==bb component + * + *

return: + * - 0 + * + * \param ndim number of dimensions of the integral (2 here) + * \param x point where the function should be evaluated + * \param ncomp number of components of the integrand (1 here) + * \param f function value + * \param userdata additional user parameters (required by the interface, NULL here) + */ +int TLinePWaveGapIntegralCuhre::Integrand_aa(const int *ndim, const double x[], + const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc} +{ + double z = x[1]*fPar[3]; + double deltasq(pow(z*fPar[1],2.0)); + f[0] = (1.0-z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + return 0; +} + +//----------------------------------------------------------------------------- +/** + *

Calculate the function value for the use with Cuhre---actual implementation of the function + * for p-wave line, cc component + * + *

return: + * - 0 + * + * \param ndim number of dimensions of the integral (2 here) + * \param x point where the function should be evaluated + * \param ncomp number of components of the integrand (1 here) + * \param f function value + * \param userdata additional user parameters (required by the interface, NULL here) + */ +int TLinePWaveGapIntegralCuhre::Integrand_cc(const int *ndim, const double x[], + const int *ncomp, double f[], void *userdata) // x = {E, z}, fPar = {twokBT, Delta(T), Ec, zc} +{ + double z = x[1]*fPar[3]; + double deltasq(pow(z*fPar[1],2.0)); + f[0] = (z*z)/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0); + return 0; +} + +//----------------------------------------------------------------------------- + std::vector TDWaveGapIntegralCuhre::fPar; /** diff --git a/src/external/BMWtools/BMWIntegrator.h b/src/external/BMWtools/BMWIntegrator.h index a901004a..3aa428c4 100644 --- a/src/external/BMWtools/BMWIntegrator.h +++ b/src/external/BMWtools/BMWIntegrator.h @@ -249,6 +249,47 @@ inline double TMCIntegrator::IntegrateFunc(size_t dim, double *x1, double *x2) return fMCIntegrator->Integral(fFunc, dim, x1, x2, (this)); } +//----------------------------------------------------------------------------- +/** + *

Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model + * assuming a cylindrical Fermi surface and a point p symmetry of the superconducting order parameter. + * The integration uses the Cuhre algorithm of the Cuba library. + */ +class TPointPWaveGapIntegralCuhre { + public: + TPointPWaveGapIntegralCuhre() : fNDim(2) {} + ~TPointPWaveGapIntegralCuhre() { fPar.clear(); } + void SetParameters(const std::vector &par) { fPar=par; } + static int Integrand_aa(const int*, const double[], const int*, double[], void*); + static int Integrand_cc(const int*, const double[], const int*, double[], void*); + double IntegrateFunc(int tag); + + protected: + static std::vector fPar; ///< parameters of the integrand + unsigned int fNDim; ///< dimension of the integral +}; + +//----------------------------------------------------------------------------- +/** + *

Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model + * assuming a cylindrical Fermi surface and a line p symmetry of the superconducting order parameter. + * The integration uses the Cuhre algorithm of the Cuba library. + */ +class TLinePWaveGapIntegralCuhre { + public: + TLinePWaveGapIntegralCuhre() : fNDim(2) {} + ~TLinePWaveGapIntegralCuhre() { fPar.clear(); } + void SetParameters(const std::vector &par) { fPar=par; } + static int Integrand_aa(const int*, const double[], const int*, double[], void*); + static int Integrand_cc(const int*, const double[], const int*, double[], void*); + double IntegrateFunc(int tag); + + protected: + static std::vector fPar; ///< parameters of the integrand + unsigned int fNDim; ///< dimension of the integral +}; + +//----------------------------------------------------------------------------- /** *

Two-dimensional integrator class for the efficient calculation of the superfluid density within the semi-classical model * assuming a cylindrical Fermi surface and a d_{x^2-y^2} symmetry of the superconducting order parameter. diff --git a/src/external/libGapIntegrals/GapIntegrals.pdf b/src/external/libGapIntegrals/GapIntegrals.pdf index b03edf6c..37334a47 100644 Binary files a/src/external/libGapIntegrals/GapIntegrals.pdf and b/src/external/libGapIntegrals/GapIntegrals.pdf differ diff --git a/src/external/libGapIntegrals/GapIntegrals.tex b/src/external/libGapIntegrals/GapIntegrals.tex index 75812971..c319c67d 100644 --- a/src/external/libGapIntegrals/GapIntegrals.tex +++ b/src/external/libGapIntegrals/GapIntegrals.tex @@ -74,16 +74,17 @@ E-Mail: & \verb?andreas.suter@psi.ch? && \section*{\musrfithead plug-in for the calculation of the temperature dependence of $\bm{1/\lambda^2}$ for various gap symmetries}% This memo is intended to give a short summary of the background on which the \gapint plug-in for \musrfit \cite{musrfit} has been developed. The aim of this implementation is the efficient calculation of integrals of the form -\begin{equation}\label{int} +\begin{equation}\label{eq:int_phi} I(T) = 1 + \frac{1}{\pi}\int_0^{2\pi}\int_{\Delta(\varphi,T)}^{\infty}\left(\frac{\partial f}{\partial E}\right) \frac{E}{\sqrt{E^2-\Delta^2(\varphi,T)}}\mathrm{d}E\mathrm{d}\varphi\,, \end{equation} where $f = (1+\exp(E/k_{\mathrm B}T))^{-1}$, like they appear e.g. in the theoretical temperature dependence of $1/\lambda^2$~\cite{Manzano}. -In order not to do too many unnecessary function calls during the final numerical evaluation we simplify the integral (\ref{int}) as far as possible analytically. The derivative of $f$ is given by +For gap symmetries which involve not only a $E$- and $\varphi$-dependence but also a $\theta$-dependence, see the special section towards the end of the memo. +In order not to do too many unnecessary function calls during the final numerical evaluation we simplify the integral (\ref{eq:int_phi}) as far as possible analytically. The derivative of $f$ is given by \begin{equation}\label{derivative} \frac{\partial f}{\partial E} = -\frac{1}{k_{\mathrm B}T}\frac{\exp(E/k_{\mathrm B}T)}{\left(1+\exp(E/k_{\mathrm B}T)\right)^2} = -\frac{1}{4k_{\mathrm B}T} \frac{1}{\cosh^2\left(E/2k_{\mathrm B}T\right)}. \end{equation} -Using (\ref{derivative}) and doing the substitution $E'^2 = E^2-\Delta^2(\varphi,T)$, equation (\ref{int}) can be written as -\begin{equation} +Using (\ref{derivative}) and doing the substitution $E'^2 = E^2-\Delta^2(\varphi,T)$, equation (\ref{eq:int_phi}) can be written as +\begin{equation}\label{eq:bmw_2d} \begin{split} I(T) & = 1 - \frac{1}{4\pi k_{\mathrm B}T}\int_0^{2\pi}\int_{\Delta(\varphi,T)}^{\infty}\frac{1}{\cosh^2\left(E/2k_{\mathrm B}T\right)}\frac{E}{\sqrt{E^2-\Delta^2(\varphi,T)}}\mathrm{d}E\mathrm{d}\varphi \\ & = 1 - \frac{1}{4\pi k_{\mathrm B}T}\int_0^{2\pi}\int_{0}^{\infty}\frac{1}{\cosh^2\left(\sqrt{E'^2+\Delta^2(\varphi,T)}/2k_{\mathrm B}T\right)}\mathrm{d}E'\mathrm{d}\varphi\,. @@ -167,6 +168,22 @@ The \gapint plug-in calculates $\tilde{I}(T)$ for the following $\Delta(\varphi) Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, $a~(1)$, $[c_0~(1),~a_{\rm G}~(1)]$. If $c_0$ and $a_{\rm G}$ are provided, the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. +\item[\textit{p}-wave (point) \cite{Pang2015}:] + \begin{equation} + \Delta(\theta, T) = \Delta(T) \sin(\theta) = \Delta(T) \cdot \sqrt{1-z^2} + \end{equation} + \musrfit theory line: \verb?userFcn libGapIntegrals TGapPointPWave 1 2 [3 [4 5]]?\\[1.5ex] + Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, [ \verb?orientation_tag?, $[c_0~(1),~a_{\rm G}~(1)]$]. If $c_0$ and $a_{\rm G}$ are provided, + the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. \\ + \verb?orientation_tag?: $0=\{aa,bb\}$, $1=cc$, and the default $2=$ average (see Eq.\ (\ref{eq:n_avg})) +\item[\textit{p}-wave (line) \cite{Ozaki1986}:] + \begin{equation} + \Delta(\theta, T) = \Delta(T) \cos(\theta) = \Delta(T) \cdot z + \end{equation} + \musrfit theory line: \verb?userFcn libGapIntegrals TGapLinePWave 1 2 [3 [4 5]]?\\[1.5ex] + Parameters: $T_{\mathrm c}~(\mathrm{K})$, $\Delta_0~(\mathrm{meV})$, [ \verb?orientation_tag?, $[c_0~(1),~a_{\rm G}~(1)]$]. If $c_0$ and $a_{\rm G}$ are provided, + the temperature dependence according to Eq.(\ref{eq:gapT_Prozorov}) will be used, otherwise Eq.(\ref{eq:gapT_Manzano}) will be utilized. \\ + \verb?orientation_tag?: $0=\{aa,bb\}$, $1=cc$, and the default $2=$ average (see Eq.\ (\ref{eq:n_avg})) \end{description} \noindent It is also possible to calculate a power law temperature dependence (in the two fluid approximation $n=4$) and the dirty \textit{s}-wave expression. @@ -193,6 +210,69 @@ Obviously for this no integration is needed. within the semi-classical model assuming a cylindrical Fermi surface and a mixed $d_{x^2-y^2} + s$ symmetry of the superconducting order parameter (effectively: $d_{x^2-y^2}$ with shifted nodes and \textit{a}-\textit{b}-anisotropy)) see the source code. +\subsection*{Gap Integrals for $\bm{\theta}$-, and $\bm{(\theta, \varphi)}$-dependent Gaps}% + +First some general formulae as found in Ref.\,\cite{Prozorov}. It assumes an anisotropic response which can be classified in 3 directions ($a$, $b$, and $c$). + +\noindent For the case of a 2D Fermi surface (cylindrical symmetry): + +\begin{equation}\label{eq:n_anisotrope_2D} +n_{aa \atop bb}(T) = 1 - \frac{1}{2\pi k_{\rm B} T} \int_0^{2\pi} \mathrm{d}\varphi\, {\cos^2(\varphi) \atop \sin^2(\varphi)} \underbrace{\int_0^\infty \mathrm{d}\varepsilon\, \left\{ \cosh\left[\frac{\sqrt{\varepsilon^2 + \Delta^2}}{2 k_{\rm B}T}\right]\right\}^{-2}}_{= G(\Delta(\varphi), T)} +\end{equation} + +\noindent For the case of a 3D Fermi surface: + +\begin{eqnarray} + n_{aa \atop bb}(T) &=& 1 - \frac{3}{4\pi k_{\rm B} T} \int_0^1 \mathrm{d}z\, (1-z^2) \int_0^{2\pi} \mathrm{d}\varphi\, {\cos^2(\varphi) \atop \sin^2(\varphi)} \cdot G(\Delta(z,\varphi), T) \label{eq:n_anisotrope_3D_aabb} \\ + n_{cc}(T) &=& 1 - \frac{3}{2\pi k_{\rm B} T} \int_0^1 \mathrm{d}z\, z^2 \int_0^{2\pi} \mathrm{d}\varphi\, \cos^2(\varphi) \cdot G(\Delta(z,\varphi), T) \label{eq:n_anisotrope_3D_cc} +\end{eqnarray} + +\noindent The ``powder averaged'' superfluid density is then defined as + +\begin{equation}\label{eq:n_avg} + n_{\rm S} = \frac{1}{3}\cdot \left[ \sqrt{n_{aa} n_{bb}} + \sqrt{n_{aa} n_{cc}} + \sqrt{n_{bb} n_{cc}} \right] +\end{equation} + +\subsubsection*{Isotropic s-Wave Gap} + +\noindent For the 2D/3D case this means that $\Delta$ is just a constant. + +\noindent For the 2D case it follows + +\begin{equation} + n_{aa \atop bb}(T) = 1 - \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T) = n_{\rm S}(T). +\end{equation} + +\noindent This is the same as Eq.(\ref{eq:bmw_2d}), assuming a $\Delta \neq f(\varphi)$. + +\vspace{5mm} + +\noindent The 3D case for $\Delta \neq f(\theta, \varphi)$: + +\noindent The variable transformation $z = \cos(\theta)$ leads to $\mathrm{d}z = -\sin(\theta)\,\mathrm{d}\theta$, $z=0 \to \theta=\pi/2$, $z=1 \to \theta=0$, and hence to + +\begin{eqnarray*} + n_{aa \atop bb}(T) &=& 1 - \frac{3}{4\pi k_{\rm B} T} \underbrace{\int_0^{\pi/2} \mathrm{d}\theta \, \sin^3(\theta)}_{= 2/3} \, \underbrace{\int_0^{2\pi} \mathrm{d}\varphi {\cos^2(\varphi) \atop \sin^2(\varphi)}}_{=\pi} \cdot G(\Delta, T) \\ + &=& 1 - \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T). \\ + n_{cc}(T) &=& 1 - \frac{3}{2\pi k_{\rm B} T} \underbrace{\int_0^{\pi/2} \mathrm{d}\theta \, \cos^2(\theta)\sin(\theta)}_{=1/3} \, \underbrace{\int_0^{2\pi} \mathrm{d}\varphi \cos^2(\varphi)}_{=\pi} \cdot G(\Delta, T) \\ + &=& 1 - \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T). +\end{eqnarray*} + +\noindent And hence + +$$ n_{\rm S}(T) = 1- \frac{1}{2 k_{\rm B} T} \cdot G(\Delta, T). $$ + +\subsubsection*{3D Fermi Surface Gap $\mathbf{\Delta \neq f(\bm\varphi)}$} + +For this case the superfluid density integrals reduce to ($z=\cos(\theta)$) + +\begin{eqnarray} + n_{aa \atop bb}(T) &=& 1 - \frac{3}{4 k_{\rm B} T} \int_0^1 \mathrm{d}z\, (1-z^2) \cdot G(\Delta(z, T),T) \\ + n_{cc}(T) &=& 1 - \frac{3}{2 k_{\rm B} T} \int_0^1 \mathrm{d}z\, z^2 \cdot G(\Delta(z, T),T) +\end{eqnarray} + + + \subsection*{License} The \gapint library is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation \cite{GPL}; either version 2 of the License, or (at your option) any later version. @@ -208,6 +288,8 @@ The \gapint library is free software; you can redistribute it and/or modify it u \bibitem{Matsui} H.~Matsui~\textit{et al.}, \textit{Direct Observation of a Nonmonotonic $d_{x^2-y^2}$-Wave Superconducting Gap in the Electron-Doped High-T$_{\mathrm c}$ Superconductor Pr$_{0.89}$LaCe$_{0.11}$CuO$_4$}, Phys.~Rev.~Lett.~\textbf{95}~(2005)~017003 \bibitem{Eremin} I.~Eremin, E.~Tsoncheva, and A.V.~Chubukov, \textit{Signature of the nonmonotonic $d$-wave gap in electron-doped cuprates}, Phys.~Rev.~B~\textbf{77}~(2008)~024508 \bibitem{AnisotropicSWave} ?? +\bibitem{Pang2015} G.M.~Pang, \emph{et al.}, Phys.~Rev.~B~\textbf{91}~(2015)~220502(R), and references in there. +\bibitem{Ozaki1986} M.~Ozaki, \emph{et al.}, Prog.~Theor.~Phys.~\textbf{75}~(1986)~442. \bibitem{Tinkham} M.~Tinkham, \textit{Introduction to Superconductivity} $2^{\rm nd}$ ed. (Dover Publications, New York, 2004). \bibitem{GPL} http://www.gnu.org/licenses/old-licenses/gpl-2.0.html diff --git a/src/external/libGapIntegrals/TGapIntegrals.cpp b/src/external/libGapIntegrals/TGapIntegrals.cpp index 702f1589..36c06122 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.cpp +++ b/src/external/libGapIntegrals/TGapIntegrals.cpp @@ -36,6 +36,8 @@ #define TWOPI 6.28318530717958647692 ClassImp(TGapSWave) +ClassImp(TGapPointPWave) +ClassImp(TGapLinePWave) ClassImp(TGapDWave) ClassImp(TGapCosSqDWave) ClassImp(TGapSinSqDWave) @@ -46,6 +48,8 @@ ClassImp(TGapPowerLaw) ClassImp(TGapDirtySWave) ClassImp(TLambdaSWave) +ClassImp(TLambdaPointPWave) +ClassImp(TLambdaLinePWave) ClassImp(TLambdaDWave) ClassImp(TLambdaAnSWave) ClassImp(TLambdaNonMonDWave1) @@ -53,6 +57,8 @@ ClassImp(TLambdaNonMonDWave2) ClassImp(TLambdaPowerLaw) ClassImp(TLambdaInvSWave) +ClassImp(TLambdaInvPointPWave) +ClassImp(TLambdaInvLinePWave) ClassImp(TLambdaInvDWave) ClassImp(TLambdaInvAnSWave) ClassImp(TLambdaInvNonMonDWave1) @@ -77,6 +83,38 @@ TGapSWave::TGapSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

point p wave gap integral + */ +TGapPointPWave::TGapPointPWave() { + TPointPWaveGapIntegralCuhre *gapint = new TPointPWaveGapIntegralCuhre(); + fGapIntegral = gapint; + gapint = nullptr; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +//-------------------------------------------------------------------- +/** + *

line p wave gap integral + */ +TGapLinePWave::TGapLinePWave() { + TLinePWaveGapIntegralCuhre *gapint = new TLinePWaveGapIntegralCuhre(); + fGapIntegral = gapint; + gapint = nullptr; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + //-------------------------------------------------------------------- /** *

@@ -181,6 +219,22 @@ TLambdaSWave::TLambdaSWave() { fLambdaInvSq = new TGapSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaPointPWave::TLambdaPointPWave() { + fLambdaInvSq = new TGapPointPWave(); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaLinePWave::TLambdaLinePWave() { + fLambdaInvSq = new TGapLinePWave(); +} + //-------------------------------------------------------------------- /** *

@@ -221,6 +275,22 @@ TLambdaInvSWave::TLambdaInvSWave() { fLambdaInvSq = new TGapSWave(); } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvPointPWave::TLambdaInvPointPWave() { + fLambdaInvSq = new TGapPointPWave(); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvLinePWave::TLambdaInvLinePWave() { + fLambdaInvSq = new TGapLinePWave(); +} + //-------------------------------------------------------------------- /** *

@@ -268,6 +338,36 @@ TGapSWave::~TGapSWave() { fPar.clear(); } +//-------------------------------------------------------------------- +/** + *

+ */ +TGapPointPWave::~TGapPointPWave() { + delete fGapIntegral; + fGapIntegral = nullptr; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TGapLinePWave::~TGapLinePWave() { + delete fGapIntegral; + fGapIntegral = nullptr; + + fTemp.clear(); + fTempIter = fTemp.end(); + fIntegralValues.clear(); + fCalcNeeded.clear(); + fPar.clear(); +} + //-------------------------------------------------------------------- /** *

@@ -367,6 +467,24 @@ TLambdaSWave::~TLambdaSWave() { fLambdaInvSq = nullptr; } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaPointPWave::~TLambdaPointPWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaLinePWave::~TLambdaLinePWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + //-------------------------------------------------------------------- /** *

@@ -412,6 +530,24 @@ TLambdaInvSWave::~TLambdaInvSWave() { fLambdaInvSq = nullptr; } +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvPointPWave::~TLambdaInvPointPWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + +//-------------------------------------------------------------------- +/** + *

+ */ +TLambdaInvLinePWave::~TLambdaInvLinePWave() { + delete fLambdaInvSq; + fLambdaInvSq = nullptr; +} + //-------------------------------------------------------------------- /** *

@@ -528,6 +664,191 @@ double TGapSWave::operator()(double t, const std::vector &par) const { } +//-------------------------------------------------------------------- +/** + *

prepare the needed parameters for the integration carried out in TPointPWaveGapIntegralCuhre. + * For details see also the Memo GapIntegrals.pdf, , especially Eq.(19) and (20). + */ +double TGapPointPWave::operator()(double t, const std::vector &par) const { + + // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] orientation tag, [[3] c0 (1), [4] aG (1)]] + + assert((par.size() >= 2) && (par.size() <= 5)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205 + // 4 or 5 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67 + // and Erratum Supercond. Sci. Technol. 21 (2008) 082003 + // c0 in the original context is c0 = (pi kB Tc) / Delta0 + // orientation tag: 0=aa,bb; 1=cc; 2=(sqrt[aa bb] + sqrt[aa cc] + sqrt[bb cc])/3 (default) + if (t <= 0.0) + return 1.0; + else if (t >= par[0]) + return 0.0; + + // check if orientation tag is given + int orientation_tag(2); + if ((par.size()==3) || (par.size()==5)) + orientation_tag = static_cast(par[2]); + + bool integralParChanged(false); + + if (fPar.empty()) { // first time calling this routine + fPar = par; + integralParChanged = true; + } else { // check if Tc or Delta0 have changed + for (unsigned int i(0); i intPar; // parameters for the integral, T & Delta(T) + intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + if ((par.size() == 2) || (par.size() == 3)) { // Carrington/Manzano + intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); + } else { // Prozorov/Giannetta + intPar.push_back(par[1]*tanh(par[2]*sqrt(par[3]*(par[0]/t-1.0)))); // Delta0*tanh(c0*sqrt(aG*(Tc/T-1))) + } + intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy + intPar.push_back(1.0); // upper limit of theta-integration + + fGapIntegral->SetParameters(intPar); + if (orientation_tag == 0) // aa,bb + ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + else if (orientation_tag == 1) // cc + ds = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + else { // average + ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + ds1 = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + ds = (ds + 2.0 * sqrt(ds*ds1))/3.0; // since aa==bb the avg looks like this + } + + intPar.clear(); + + if (newTemp) + fIntegralValues.push_back(ds); + else + fIntegralValues[vectorIndex] = ds; + + fCalcNeeded[vectorIndex] = false; + } + + return fIntegralValues[vectorIndex]; +} + +//-------------------------------------------------------------------- +/** + *

prepare the needed parameters for the integration carried out in TLinePWaveGapIntegralCuhre. + * For details see also the Memo GapIntegrals.pdf, especially Eq.(19) and (20). + */ +double TGapLinePWave::operator()(double t, const std::vector &par) const { + + // parameters: [0] Tc (K), [1] Delta0 (meV), [[2] orientation tag, [[3] c0 (1), [4] aG (1)]] + + assert((par.size() >= 2) && (par.size() <= 5)); // 2 parameters: see A.~Carrington and F.~Manzano, Physica~C~\textbf{385}~(2003)~205 + // 4 or 5 parameters: see R. Prozorov and R. Giannetta, Supercond. Sci. Technol. 19 (2006) R41-R67 + // and Erratum Supercond. Sci. Technol. 21 (2008) 082003 + // c0 in the original context is c0 = (pi kB Tc) / Delta0 + // orientation tag: 0=aa,bb; 1=cc; 2=(sqrt[aa bb] + sqrt[aa cc] + sqrt[bb cc])/3 (default) + if (t <= 0.0) + return 1.0; + else if (t >= par[0]) + return 0.0; + + + // check if orientation tag is given + int orientation_tag(2); + if ((par.size()==3) || (par.size()==5)) + orientation_tag = static_cast(par[2]); + + bool integralParChanged(false); + + if (fPar.empty()) { // first time calling this routine + fPar = par; + integralParChanged = true; + } else { // check if parameter have changed + for (unsigned int i(0); i intPar; // parameters for the integral, T & Delta(T) + intPar.push_back(0.172346648*t); // 2 kB T, kB in meV/K = 0.086173324 meV/K + if ((par.size() == 2) || (par.size() == 3)) { // Carrington/Manzano + intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51))); + } else { // Prozorov/Giannetta + intPar.push_back(par[1]*tanh(par[3]*sqrt(par[4]*(par[0]/t-1.0)))); // Delta0*tanh(c0*sqrt(aG*(Tc/T-1))) + } + intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy + intPar.push_back(1.0); // upper limit of z-integration + + fGapIntegral->SetParameters(intPar); + if (orientation_tag == 0) // aa,bb + ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + else if (orientation_tag == 1) // cc + ds = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + else { // average + ds = 1.0-(intPar[2]*3.0)/(2.0*intPar[0])*fGapIntegral->IntegrateFunc(0); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + ds1 = 1.0-(intPar[2]*3.0)/(intPar[0])*fGapIntegral->IntegrateFunc(1); // integral prefactor is by 2 lower [Eqs.(19,20)] since intPar[0]==2kB T! + ds = (ds + 2.0 * sqrt(ds*ds1))/3.0; // since aa==bb the avg looks like this + } + + intPar.clear(); + + if (newTemp) + fIntegralValues.push_back(ds); + else + fIntegralValues[vectorIndex] = ds; + + fCalcNeeded[vectorIndex] = false; + } + + return fIntegralValues[vectorIndex]; +} + //-------------------------------------------------------------------- /** *

prepare the needed parameters for the integration carried out in TDWaveGapIntegralCuhre. @@ -609,7 +930,6 @@ double TGapDWave::operator()(double t, const std::vector &par) const { } return fIntegralValues[vectorIndex]; - } //-------------------------------------------------------------------- @@ -1103,7 +1423,40 @@ double TLambdaSWave::operator()(double t, const std::vector &par) const return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); +} +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaPointPWave::operator()(double t, const std::vector &par) const +{ + assert(par.size() == 2); // two parameters: Tc, Delta0 + + if (t >= par[0]) + return -1.0; + + if (t <= 0.0) + return 1.0; + + return 1.0/sqrt((*fLambdaInvSq)(t, par)); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaLinePWave::operator()(double t, const std::vector &par) const +{ + assert(par.size() == 2); // two parameters: Tc, Delta0 + + if (t >= par[0]) + return -1.0; + + if (t <= 0.0) + return 1.0; + + return 1.0/sqrt((*fLambdaInvSq)(t, par)); } //-------------------------------------------------------------------- @@ -1121,7 +1474,6 @@ double TLambdaDWave::operator()(double t, const std::vector &par) const return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1139,7 +1491,6 @@ double TLambdaAnSWave::operator()(double t, const std::vector &par) cons return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1157,7 +1508,6 @@ double TLambdaNonMonDWave1::operator()(double t, const std::vector &par) return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1175,7 +1525,6 @@ double TLambdaNonMonDWave2::operator()(double t, const std::vector &par) return 1.0; return 1.0/sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1192,7 +1541,6 @@ double TLambdaPowerLaw::operator()(double t, const std::vector &par) con return -1.0; return 1.0/sqrt(1.0 - pow(t/par[0], par[1])); - } @@ -1211,7 +1559,40 @@ double TLambdaInvSWave::operator()(double t, const std::vector &par) con return 1.0; return sqrt((*fLambdaInvSq)(t, par)); +} +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaInvPointPWave::operator()(double t, const std::vector &par) const +{ + assert(par.size() == 2); // two parameters: Tc, Delta0 + + if (t >= par[0]) + return 0.0; + + if (t <= 0.0) + return 1.0; + + return sqrt((*fLambdaInvSq)(t, par)); +} + +//-------------------------------------------------------------------- +/** + *

+ */ +double TLambdaInvLinePWave::operator()(double t, const std::vector &par) const +{ + assert(par.size() == 2); // two parameters: Tc, Delta0 + + if (t >= par[0]) + return 0.0; + + if (t <= 0.0) + return 1.0; + + return sqrt((*fLambdaInvSq)(t, par)); } //-------------------------------------------------------------------- @@ -1229,7 +1610,6 @@ double TLambdaInvDWave::operator()(double t, const std::vector &par) con return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1247,7 +1627,6 @@ double TLambdaInvAnSWave::operator()(double t, const std::vector &par) c return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1265,7 +1644,6 @@ double TLambdaInvNonMonDWave1::operator()(double t, const std::vector &p return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1283,7 +1661,6 @@ double TLambdaInvNonMonDWave2::operator()(double t, const std::vector &p return 1.0; return sqrt((*fLambdaInvSq)(t, par)); - } //-------------------------------------------------------------------- @@ -1300,7 +1677,6 @@ double TLambdaInvPowerLaw::operator()(double t, const std::vector &par) return 0.0; return sqrt(1.0 - pow(t/par[0], par[1])); - } //-------------------------------------------------------------------- diff --git a/src/external/libGapIntegrals/TGapIntegrals.h b/src/external/libGapIntegrals/TGapIntegrals.h index a0d344ee..07c498be 100644 --- a/src/external/libGapIntegrals/TGapIntegrals.h +++ b/src/external/libGapIntegrals/TGapIntegrals.h @@ -63,6 +63,62 @@ private: ClassDef(TGapSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ +class TGapPointPWave : public PUserFcnBase { + +public: + TGapPointPWave(); + virtual ~TGapPointPWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TPointPWaveGapIntegralCuhre *fGapIntegral; + mutable std::vector fTemp; + mutable std::vector::const_iterator fTempIter; + mutable std::vector fIntegralValues; + mutable std::vector fCalcNeeded; + + mutable std::vector fPar; + + ClassDef(TGapPointPWave,1) +}; + +//-------------------------------------------------------------------- +/** + *

+ */ +class TGapLinePWave : public PUserFcnBase { + +public: + TGapLinePWave(); + virtual ~TGapLinePWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TLinePWaveGapIntegralCuhre *fGapIntegral; + mutable std::vector fTemp; + mutable std::vector::const_iterator fTempIter; + mutable std::vector fIntegralValues; + mutable std::vector fCalcNeeded; + + mutable std::vector fPar; + + ClassDef(TGapLinePWave,1) +}; + //-------------------------------------------------------------------- /** *

@@ -297,6 +353,50 @@ private: ClassDef(TLambdaSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaPointPWave : public PUserFcnBase { + +public: + TLambdaPointPWave(); + virtual ~TLambdaPointPWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapPointPWave *fLambdaInvSq; + + ClassDef(TLambdaPointPWave,1) +}; + +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaLinePWave : public PUserFcnBase { + +public: + TLambdaLinePWave(); + virtual ~TLambdaLinePWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapLinePWave *fLambdaInvSq; + + ClassDef(TLambdaLinePWave,1) +}; + //-------------------------------------------------------------------- /** *

@@ -428,6 +528,50 @@ private: ClassDef(TLambdaInvSWave,1) }; +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaInvPointPWave : public PUserFcnBase { + +public: + TLambdaInvPointPWave(); + virtual ~TLambdaInvPointPWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapPointPWave *fLambdaInvSq; + + ClassDef(TLambdaInvPointPWave,1) +}; + +//-------------------------------------------------------------------- +/** + *

+ */ +class TLambdaInvLinePWave : public PUserFcnBase { + +public: + TLambdaInvLinePWave(); + virtual ~TLambdaInvLinePWave(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(std::vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const std::vector&) const; + +private: + TGapLinePWave *fLambdaInvSq; + + ClassDef(TLambdaInvLinePWave,1) +}; + //-------------------------------------------------------------------- /** *

diff --git a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h index 66ee9a3d..37e17f68 100644 --- a/src/external/libGapIntegrals/TGapIntegralsLinkDef.h +++ b/src/external/libGapIntegrals/TGapIntegralsLinkDef.h @@ -35,6 +35,8 @@ #pragma link off all functions; #pragma link C++ class TGapSWave+; +#pragma link C++ class TGapPointPWave+; +#pragma link C++ class TGapLinePWave+; #pragma link C++ class TGapDWave+; #pragma link C++ class TGapCosSqDWave+; #pragma link C++ class TGapSinSqDWave+; @@ -44,12 +46,16 @@ #pragma link C++ class TGapPowerLaw+; #pragma link C++ class TGapDirtySWave+; #pragma link C++ class TLambdaSWave+; +#pragma link C++ class TLambdaPointPWave+; +#pragma link C++ class TLambdaLinePWave+; #pragma link C++ class TLambdaDWave+; #pragma link C++ class TLambdaAnSWave+; #pragma link C++ class TLambdaNonMonDWave1+; #pragma link C++ class TLambdaNonMonDWave2+; #pragma link C++ class TLambdaPowerLaw+; #pragma link C++ class TLambdaInvSWave+; +#pragma link C++ class TLambdaInvPointPWave+; +#pragma link C++ class TLambdaInvLinePWave+; #pragma link C++ class TLambdaInvDWave+; #pragma link C++ class TLambdaInvAnSWave+; #pragma link C++ class TLambdaInvNonMonDWave1+;