Merge pull request #127 from ytdHuang/yide-patch-1

Make `BosonBath` align with `QuTiP-BoFiN`

.github/ISSUE_TEMPLATE/bug_report.yaml

@@ -0,0 +1,53 @@
+name: 🐛 Bug Report
+description: Spotted a bug? Report it to us!
+labels: ["bug"]
+body:
+  - type: textarea
+    id: bug-description
+    attributes:
+      label: Bug Description
+      description: Tell us what went wrong (including what triggered the bug)
+      placeholder: "A clear and concise description of what the bug is, and the steps to reproduce it"
+    validations:
+      required: true
+  - type: textarea
+    id: code-to-reproduce
+    attributes:
+      label: Code to Reproduce the Bug
+      description: Please provide a minimal working example. Paste your code directly (It will be automatically formatted, so there's no need for backticks)
+      placeholder: "using HierarchicalEOM\nM_S(sigmax())"
+      render: julia
+  - type: textarea
+    id: bug-output
+    attributes:
+      label: Code Output
+      description: Please paste the relevant output here (automatically formatted)
+      placeholder: "Schrodinger Eq. type HEOMLS matrix acting on even-parity ADOs\nsystem dims = [2]\nnumber of ADOs N = 1\ndata =\nMatrixOperator(4 × 4)"
+      render: shell
+  - type: textarea
+    id: expected-behaviour
+    attributes:
+      label: Expected Behaviour
+      description: What should have happened instead?
+      placeholder: "A clear and concise description of what you expected to happen"
+    validations:
+      required: true
+  - type: textarea
+    id: your-environment
+    attributes:
+      label: Your Environment
+      description: Please use `HierarchicalEOM.about()` or `HierarchicalEOM.versioninfo()` to get the information about your environment and paste it here (automatically formatted)
+      placeholder: "Julia              Ver. ***\nHierarchicalEOM    Ver. ***\nQuantumToolbox     Ver. ***\nSciMLOperators     Ver. ***\nLinearSolve        Ver. ***\nOrdinaryDiffEqCore Ver. ***\nOS       : ***\nWORD_SIZE: ***\nLIBM     : ***\nLLVM     : ***\nBLAS     : ***"
+      render: shell
+    validations:
+      required: true
+  - type: textarea
+    id: additional-context
+    attributes:
+      label: Additional Context
+      description: Add anything else you want to tell us.
+      placeholder: "You can include output of `using Pkg; Pkg.status()`, `using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)`, or even images, plots, etc."
+  - type: markdown
+    id: thank-you
+    attributes:
+      value: Thanks for taking the time to fill out this bug report!

.github/ISSUE_TEMPLATE/feature_request.yaml

@@ -0,0 +1,36 @@
+name: 🚀 Feature Request
+description: Suggest an idea for this project!
+labels: ["enhancement"]
+body:
+  - type: textarea
+    id: problem-description
+    attributes:
+      label: Problem Description
+      description: Give us a clear and concise description of what the problem is!
+      placeholder: "Example - It would be better if [...]"
+    validations:
+      required: true
+  - type: textarea
+    id: proposed-solution
+    attributes:
+      label: Proposed Solution
+      description: Give us a description of your proposed solution
+      placeholder: "A clear and concise description of what you want to happen"
+    validations:
+      required: true
+  - type: textarea
+    id: alternate-solutions
+    attributes:
+      label: Alternate Solutions
+      description: Are there other ways of implementing this feature?
+      placeholder: "A clear and concise description of any alternative solutions or features you've considered"
+  - type: textarea
+    id: additional-context
+    attributes:
+      label: Additional Context
+      description: You can tell us more about your idea, add code snippets or images, etc.
+      placeholder: "Give us more context or screenshots about the feature request here"
+  - type: markdown
+    id: thank-you
+    attributes:
+      value: Thanks for taking the time to fill out this feature request!

.github/ISSUE_TEMPLATE/others.yaml

@@ -0,0 +1,15 @@
+name: 🤔 Others
+description: Tell us about issues that aren't exactly bugs or features
+body:
+  - type: textarea
+    id: issue-description
+    attributes:
+      label: Describe the Issue!
+      description: Tell us about your issue. You can include screenshots, code snippets, etc.
+      placeholder: "Please describe the issue here"
+    validations:
+      required: true
+  - type: markdown
+    id: thank-you
+    attributes:
+      value: Thanks for taking the time to fill out this issue!

.github/pull_request_template.md

@@ -0,0 +1,18 @@
+## Checklist
+Thank you for contributing to `HierarchicalEOM.jl`! Please make sure you have finished the following tasks before opening the PR.
+
+- [ ] Please read [Contributing to QuantumToolbox.jl](https://qutip.org/QuantumToolbox.jl/stable/resources/contributing).
+- [ ] Any code changes were done in a way that does not break public API.
+- [ ] Appropriate tests were added and tested locally by running: `make test`.
+- [ ] Any code changes should be `julia` formatted by running: `make format`.
+- [ ] All documents (in `docs/` folder) related to code changes were updated and able to build locally by running: `make docs`.
+
+Request for a review after you have completed all the tasks. If you have not finished them all, you can also open a [Draft Pull Request](https://github.blog/2019-02-14-introducing-draft-pull-requests/) to let the others know this on-going work.
+
+## Description
+Describe the proposed change here.
+
+## Related issues or PRs
+Please mention the related issues or PRs here. If the PR fixes an issue, use the keyword close/closes/closed/fix/fixes/fixed/resolve/resolves/resolved followed by the issue id, e.g. fix #[id]
+
+## Additional context

Project.toml

@@ -1,7 +1,7 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 authors = ["Yi-Te Huang <[email protected]>"]
-version = "2.3.0"
+version = "2.3.1"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
@@ -37,7 +37,7 @@ LinearSolve = "2.4.2 - 2"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1"
 Pkg = "1"
-QuantumToolbox = "0.22"
+QuantumToolbox = "0.22 - 0.23"
 Reexport = "1"
 SciMLBase = "2"
 SciMLOperators = "0.3"

README.md

@@ -90,3 +90,9 @@ If you like `HierarchicalEOM.jl`, we would appreciate it if you starred the repo
 
 ## License
 `HierarchicalEOM.jl` is released under the [BSD 3-Clause License](./LICENSE.md).
+
+## Contributing to HierarchicalEOM.jl
+
+You are most welcome to contribute to `HierarchicalEOM.jl` development by forking this repository and sending pull requests (PRs), or filing bug reports at the issues page. You can also help out with users' questions, or discuss proposed changes in the [QuTiP discussion group](https://groups.google.com/g/qutip).
+
+For more information about contribution, including technical advice, please see the [Contributing to QuantumToolbox.jl](https://qutip.org/QuantumToolbox.jl/stable/resources/contributing) section of the `QuantumToolbox.jl` documentation.

docs/src/bath_boson/Boson_Drude_Lorentz.md

@@ -1,6 +1,6 @@
 # [Drude-Lorentz Spectral Density](@id Boson-Drude-Lorentz)
 ```math
-J(\omega)=\frac{4\Delta W\omega}{\omega^2+W^2}
+J(\omega) = 2 \pi \sum_k |g_k|^2 \delta(\omega-\omega_k) = \frac{4\Delta W\omega}{\omega^2+W^2}
 ```
 Here, ``\Delta`` represents the coupling strength between system and the bosonic environment with band-width ``W``.
 

docs/src/bath_boson/Boson_Underdamped.md

@@ -1,8 +1,8 @@
 # [Underdamped Spectral Density](@id Boson-Underdamped)
 ```math
-J(\omega)=\frac{\lambda^2 \Gamma \omega}{(\omega^2 - \omega_0^2)^2 + \omega^2\Gamma^2}
+J(\omega) = 2 \pi \sum_k |g_k|^2 \delta(\omega-\omega_k) = \frac{2 \Delta^2 W \omega}{(\omega^2 - \omega_0^2)^2 + \omega^2 W^2}
 ```
-Here, ``\lambda`` represents the coupling strength between system and the bosonic environment with band-width ``\Gamma`` and resonance frequency ``\omega_0``.
+Here, ``\Delta`` represents the coupling strength between system and the bosonic environment with band-width ``W`` and resonance frequency ``\omega_0``.
 
 ## Matsubara Expansion
 With Matsubara Expansion, the correlation function can be analytically solved and expressed as follows:
@@ -12,26 +12,26 @@ C(t_1, t_2) = C^\mathrm{R}(t_1, t_2) + iC^\mathrm{I}(t_1, t_2) = \sum_{l=1}^{\in
 with
 ```math
 \begin{aligned}
-\gamma_{1}^\mathrm{R} &= -i\Omega + \frac{\Gamma}{2},\\
-\eta_{1}^\mathrm{R} &= \frac{\lambda^2}{4\Omega}\coth\left[\frac{1}{2 k_B T}\left(\Omega + i\frac{\Gamma}{2}\right)\right],\\
-\gamma_{2}^\mathrm{R} &= i\Omega + \frac{\Gamma}{2},\\
-\eta_{2}^\mathrm{R} &= \frac{\lambda^2}{4\Omega}\coth\left[\frac{1}{2 k_B T}\left(\Omega - i\frac{\Gamma}{2}\right)\right],\\
-\gamma_{l\neq 2}^\mathrm{R} &= 2\pi l k_B T,\\
-\eta_{l\neq 2}^\mathrm{R} &= -2 k_B T \cdot \frac{\lambda^2 \Gamma \cdot \gamma_l^\mathrm{R}}{\left[\left(\Omega + i\frac{\Gamma}{2}\right)^2 + {\gamma_l^\mathrm{R}}^2\right]\left[\left(\Omega - i\frac{\Gamma}{2}\right)^2 + {\gamma_l^\mathrm{R}}^2\right]},\\
-\gamma_{1}^\mathrm{I} &= i\Omega + \frac{\Gamma}{2},\\
-\eta_{1}^\mathrm{I} &= i\frac{\lambda^2}{4\Omega},\\
-\gamma_{2}^\mathrm{I} &= -i\Omega + \frac{\Gamma}{2},\\
-\eta_{2}^\mathrm{I} &= -i\frac{\lambda^2}{4\Omega},
+\gamma_{1}^\mathrm{R} &= -i\Omega + \frac{W}{2},\\
+\eta_{1}^\mathrm{R} &= \frac{\Delta^2}{4\Omega}\coth\left[\frac{1}{2 k_B T}\left(\Omega + i\frac{W}{2}\right)\right],\\
+\gamma_{2}^\mathrm{R} &= i\Omega + \frac{W}{2},\\
+\eta_{2}^\mathrm{R} &= \frac{\Delta^2}{4\Omega}\coth\left[\frac{1}{2 k_B T}\left(\Omega - i\frac{W}{2}\right)\right],\\
+\gamma_{l}^\mathrm{R} &= 2\pi l k_B T ~~\forall~~ l \geq 3,\\
+\eta_{l}^\mathrm{R} &= -2 k_B T \cdot \frac{\Delta^2 W \cdot \gamma_l^\mathrm{R}}{\left[\left(\Omega + i\frac{W}{2}\right)^2 + {\gamma_l^\mathrm{R}}^2\right]\left[\left(\Omega - i\frac{W}{2}\right)^2 + {\gamma_l^\mathrm{R}}^2\right]} ~~\forall~~ l \geq 3,\\
+\gamma_{1}^\mathrm{I} &= i\Omega + \frac{W}{2},\\
+\eta_{1}^\mathrm{I} &= i\frac{\Delta^2}{4\Omega},\\
+\gamma_{2}^\mathrm{I} &= -i\Omega + \frac{W}{2},\\
+\eta_{2}^\mathrm{I} &= -i\frac{\Delta^2}{4\Omega},
 \end{aligned}
 ```
-where ``\Omega = \sqrt{\omega_0^2 + (\Gamma/2)^2}``.
+where ``\Omega = \sqrt{\omega_0^2 - (W/2)^2}``.
 This can be constructed by the built-in function [`Boson_Underdamped_Matsubara`](@ref):
 ```julia
 Vs # coupling operator
-λ  # coupling strength
-Γ  # band-width of the environment
+Δ  # coupling strength
+W  # band-width of the environment
 ω0 # resonance frequency of the environment
 kT # the product of the Boltzmann constant k and the absolute temperature T
 N  # Number of exponential terms
-bath = Boson_Underdamped_Matsubara(Vs, λ, Γ, ω0, kT, N - 2)
+bath = Boson_Underdamped_Matsubara(Vs, Δ, W, ω0, kT, N - 2)
 ```

docs/src/bath_fermion/Fermion_Lorentz.md

@@ -1,6 +1,6 @@
 # [Lorentz Spectral Density](@id doc-Fermion-Lorentz)
 ```math
-J(\omega)=\frac{\Gamma W^2}{(\omega-\mu)^2+W^2}
+J(\omega) = 2 \pi \sum_k |g_k|^2 \delta(\omega-\omega_k) = \frac{\Gamma W^2}{(\omega-\mu)^2+W^2}
 ```
 Here, ``\Gamma`` represents the coupling strength between system and the fermionic environment with chemical potential ``\mu`` and band-width ``W``.
 

src/bath/BosonBath.jl

@@ -180,15 +180,15 @@ function BosonBath(
         deleteat!(ηI, imag_idx)
         deleteat!(γI, imag_idx)
 
+        bRI = bosonRealImag(_op, ηRI_real, ηRI_imag, γRI)
         bR = bosonReal(_op, ηR, γR)
         bI = bosonImag(_op, ηI, γI)
-        bRI = bosonRealImag(_op, ηRI_real, ηRI_imag, γRI)
-        Nterm_new = bR.Nterm + bI.Nterm + bRI.Nterm
+        Nterm_new = bRI.Nterm + bR.Nterm + bI.Nterm
         if Nterm != (Nterm_new + bRI.Nterm)
             error("Conflicts occur in combining real and imaginary parts of bath correlation function.")
         end
 
-        return BosonBath(AbstractBosonBath[bR, bI, bRI], _op, Nterm_new, δ)
+        return BosonBath(AbstractBosonBath[bRI, bR, bI], _op, Nterm_new, δ)
 
     else
         bR = bosonReal(_op, η_real, γ_real)

src/bath_correlation_functions/boson/Underdamped.jl

@@ -1,32 +1,35 @@
 export Boson_Underdamped_Matsubara
 
 @doc raw"""
-    Boson_Underdamped_Matsubara(op, λ, Γ, ω0, kT, N)
+    Boson_Underdamped_Matsubara(op, λ, W, ω0, kT, N)
 Construct an underdamped bosonic bath with Matsubara expansion
 
 # Parameters
 - `op` : The system coupling operator, must be Hermitian and, for fermionic systems, even-parity to be compatible with charge conservation.
 - `λ::Real`: The coupling strength between the system and the bath.
-- `Γ::Real`: The band-width of the bath spectral density.
+- `W::Real`: The band-width of the bath spectral density.
 - `ω0::Real`: The resonance frequency of the bath spectral density.
 - `kT::Real`: The product of the Boltzmann constant ``k`` and the absolute temperature ``T`` of the bath.
 - `N::Int`: (N+2)-terms of exponential terms are used to approximate the bath correlation function.
 
 # Returns
 - `bath::BosonBath` : a bosonic bath object with describes the interaction between system and bosonic bath
 """
-function Boson_Underdamped_Matsubara(op, λ::Real, Γ::Real, ω0::Real, kT::Real, N::Int)
-    Ω = sqrt(ω0^2 - (Γ / 2)^2)
-    ν = 2π * kT * (1:N)
+function Boson_Underdamped_Matsubara(op, λ::Real, W::Real, ω0::Real, kT::Real, N::Int)
+    Ω = sqrt(ω0^2 - (W / 2)^2)
+    ν = (2 * π * kT) .* (1:N)
 
-    η_real = ComplexF64[(λ^2/(4*Ω))*coth((Ω + im * Γ / 2) / (2 * kT)), (λ^2/(4*Ω))*coth((Ω - im * Γ / 2) / (2 * kT))]
-    γ_real = ComplexF64[Γ/2-im*Ω, Γ/2+im*Ω]
-    η_imag = ComplexF64[(λ^2/(4*Ω))*im, -(λ^2 / (4 * Ω))*im]
-    γ_imag = ComplexF64[Γ/2-im*Ω, Γ/2+im*Ω]
+    η_real = ComplexF64[(λ^2/(4*Ω))*coth((Ω + 1im * W / 2) / (2 * kT)), (λ^2/(4*Ω))*coth((Ω - 1im * W / 2) / (2 * kT))]
+    γ_real = ComplexF64[W/2-1im*Ω, W/2+1im*Ω]
+    η_imag = ComplexF64[(λ^2/(4*Ω))*1im, -1im*(λ^2/(4*Ω))]
+    γ_imag = ComplexF64[W/2-1im*Ω, W/2+1im*Ω]
 
     if N > 0
         for l in 1:N
-            append!(η_real, -2 * λ^2 * Γ * kT * ν[l] / (((Ω + im * Γ / 2)^2 + ν[l]^2) * ((Ω - im * Γ / 2)^2 + ν[l]^2)))
+            append!(
+                η_real,
+                -2 * λ^2 * W * kT * ν[l] / (((Ω + 1im * W / 2)^2 + ν[l]^2) * ((Ω - 1im * W / 2)^2 + ν[l]^2)),
+            )
             append!(γ_real, ν[l])
         end
     end

test/bath.jl

@@ -18,15 +18,17 @@
 
     ## check for combine
     b = BosonBath(op, η1, γ1)
+    types = ["bRI", "bRI", "bRI", "bRI"]
     @test length(b) == 4
 
     ## check for η and γ list, and coupling operator
-    η = []
-    γ = []
-    for e in b
-        push!(η, e.η)
-        push!(γ, e.γ)
+    η = Vector{ComplexF64}(undef, length(b))
+    γ = Vector{ComplexF64}(undef, length(b))
+    for (i, e) in enumerate(b)
+        η[i] = e.η
+        γ[i] = e.γ
         @test e.op == op
+        @test e.types == types[i]
     end
     @test η == [1.0 + 0.0im, 10.0 + 0.0im, 5.0 + 0.0im, 9.0 + 0.0im]
     @test γ == [0.1 + 0.0im, 0.3 + 0.0im, 0.5 + 0.0im, 0.7 + 0.0im]
@@ -37,20 +39,21 @@
 
     ## check for combine
     b = BosonBath(op, η1, γ1, η2, γ2)
+    types = ["bRI", "bRI", "bR", "bR", "bI", "bI", "bI"]
     @test length(b) == 7
     @test C(b, [0.183183])[1] ≈ C(bs, [0.183183])[1]
 
     ## check for η and γ list, and coupling operator
-    η = []
-    γ = []
-    for e in b
-        push!(η, e.η)
-        push!(γ, e.γ)
+    η = Vector{ComplexF64}(undef, length(b))
+    γ = Vector{ComplexF64}(undef, length(b))
+    for (i, e) in enumerate(b)
+        η[i] = e.η
+        γ[i] = e.γ
         @test e.op == op
-        @test (e.types == "bR") || (e.types == "bI") || (e.types == "bRI")
+        @test e.types == types[i]
     end
-    @test η == [10.0 + 0.0im, 9.0 + 0.0im, 4.0 + 0.0im, 8.0 + 0.0im, 10.0 + 0.0im, 1.0 + 2.0im, 5.0 + 6.0im]
-    @test γ == [0.3 + 0.0im, 0.7 + 0.0im, 0.2 + 0.0im, 0.6 + 0.0im, 0.9 + 0.0im, 0.1 + 0.0im, 0.5 + 0.0im]
+    @test η == [1.0 + 2.0im, 5.0 + 6.0im, 10.0 + 0.0im, 9.0 + 0.0im, 4.0 + 0.0im, 8.0 + 0.0im, 10.0 + 0.0im]
+    @test γ == [0.1 + 0.0im, 0.5 + 0.0im, 0.3 + 0.0im, 0.7 + 0.0im, 0.2 + 0.0im, 0.6 + 0.0im, 0.9 + 0.0im]
     @test show(devnull, MIME("text/plain"), b) === nothing
 
     ## check for exponents
@@ -75,19 +78,20 @@
     ################################################
     # Boson bath (RWA)
     b = BosonBathRWA(op, η1, γ3, η2, γ4)
+    types = ["bA", "bA", "bA", "bA", "bA", "bE", "bE", "bE", "bE", "bE"]
     @test length(bs) == 10
     cp, cm = C(b, [0.183183])
     @test cp[1] ≈ 22.384506765987076 + 0.7399082821797519im
     @test cm[1] ≈ 26.994911851482776 - 0.799138487523946im
 
     ## check for η and γ list, and coupling operator
-    η = []
-    γ = []
-    for e in b
-        push!(η, e.η)
-        push!(γ, e.γ)
+    η = Vector{ComplexF64}(undef, length(b))
+    γ = Vector{ComplexF64}(undef, length(b))
+    for (i, e) in enumerate(b)
+        η[i] = e.η
+        γ[i] = e.γ
         @test e.op == op
-        @test (e.types == "bA") || (e.types == "bE")
+        @test e.types == types[i]
     end
     @test η == [1, 3, 5, 7, 9, 2, 4, 6, 8, 10]
     @test γ == [
@@ -127,13 +131,14 @@
     ################################################
     # Fermion bath
     b = FermionBath(op, η1, γ3, η2, γ4)
+    types = ["fA", "fA", "fA", "fA", "fA", "fE", "fE", "fE", "fE", "fE"]
     @test length(b) == 10
     cp, cm = C(b, [0.183183])
     @test cp[1] ≈ 22.384506765987076 + 0.7399082821797519im
     @test cm[1] ≈ 26.994911851482776 - 0.799138487523946im
-    for e in b
+    for (i, e) in enumerate(b)
         @test e.op == op
-        @test (e.types == "fA") || (e.types == "fE")
+        @test e.types == types[i]
     end
     @test show(devnull, MIME("text/plain"), b) === nothing