Unable to understand how to fix this Julia code to have the right output

Question

I am struggling since last month with a snippet inside a code that should realize my master thesis and I couldn't find what goes wrong. But first, Just few lines for understanding it. This snippet deals with an Information measure called Quantum Conditional Mutual Information (QCMI), defined as follows for a tripartite quantum system $A \otimes B \otimes C$ :

$$

I(A:B|C) = S(\rho_{AC}) + S(\rho_{BC}) - S(\rho_{C}) - S(\rho_{ABC})

$$

As you can see, it is written as von Neumann entropy terms . Given a quantum state described by its density matrix $\rho = \ket{\psi}\bra{\psi}$,

$\psi$ being the wavefunction associated to the quantum state, the von Neumann entropy is defined as

$$

S(\rho) = \, -\, Tr(\rho\log\rho) = - \sum_{i=1}^d \lambda_{i,\rho}\,\log\lambda_{i,\rho}

$$

where the spectral decomposition of a d-dimensional quantum state has been considered , that is

$$

\rho = \sum_{i=1}^d\,\lambda_{i,\rho}\ket{\psi_i}\bra{\psi_i}

$$

and the several $\lambda_i$ are non-negative real numbers satisfying the conditions

$$

\lambda_{1,\rho} \ge \lambda_{2,\rho} \ge \dots \ge \lambda_{d,\rho}, \qquad \sum\lambda_{i,\rho} = 1

$$

The snippet is supposed to plot the QCMI $I(A:B|C)$ versus $p \in [0,1]$ for the following mixed state built on three qubits :

\begin{equation}
    \begin{split}
        \rho & = p\,\rho_{W} + (1 - p) (\rho_{EPR} \otimes \rho_0) \\
        & = p\,\ket{W}\bra{W} + (1 - p)\Bigl(\ket{\psi_{EPR}}\bra{\psi_{EPR}} \otimes \ket{0}\bra{0}\Bigr)
    \end{split}
\end{equation}

In this case,

\mathbf{\ket{W}} = \begin{pmatrix}
     0 \\
    \frac{1}{\sqrt{3}} \\
    \frac{1}{\sqrt{3}} \\
    0 \\
    \frac{1}{\sqrt{3}} \\
    0 \\
    0 \\
    0

\end{pmatrix}

has entanglement distributed among A,B and C, whereas the other term is a state for which AB is not entangled with C. It is made of two parts:

The Bell state

\mathbf{\ket{\psi_{EPR}}} = \begin{pmatrix}
    \frac{1}{\sqrt{2}} \\
    0 \\
    0 \\
    \frac{1}{\sqrt{2}}
\end{pmatrix}

and the zero-state

\mathbf{\ket{0}} = \begin{pmatrix}
    0 \\
    1
\end{pmatrix}

Here it follows what I tried to explain :

using LinearAlgebra
using Plots
using QuantumInformation
#using OpenQuantumTools

ψ0 = [0,1]

function epr_state()
    ket = zeros(4)
    ket[1] = 1/√2
    ket[4] = 1/√2
    return ket
end

function W_state()
    ket = zeros(8)
    ket[2] = 1/√3    
    ket[3] = 1/√3
    ket[5] = 1/√3
    return ket
end
function density_matrix(ket)
    return ket * ket'
end
ψ_EPR = epr_state()
ρ_EPR = density_matrix(ψ_EPR) #Density matrix |EPR><EPR|
ρ0 = density_matrix(ψ0) #Density matrix |0><0|
ψW = W_state()
ρ_W = density_matrix(ψW) #Density matrix |W><W|
#ρ = p * ρW + (1 - p)* kron(ρEPR, ρ0), where we defined the normalized density operators with unitary trace :

ρEPR = ρ_EPR/tr(ρ_EPR)

ρW = ρ_W/tr(ρ_W)

# Define the von Neumann entropy function
function von_neumann_entropy(rho)
    eigenvalues = eigvals(rho)
     entropy = sum(-eig * log2(eig) for eig in eigenvalues if eig > 0)
    #entropy = sum(-xlogx(eig)/log(2) for eig in eigenvalues if eig > 0)
    return entropy
end

function QCMI(p)
    
ψ_EPR = epr_state()
ρ_EPR = density_matrix(ψ_EPR) #Density matrix |EPR><EPR|
ρ0 = density_matrix(ψ0) #Density matrix |0><0|
ψW = W_state()
ρ_W = density_matrix(ψW) #Density matrix |W><W|
ρEPR = ρ_EPR/tr(ρ_EPR)
ρW = ρ_W/tr(ρ_W)

    ρ = p* ρW + (1 - p)* kron(ρEPR,ρ0)
    ρ_C = ptrace(ρ,[2,2,2],[1,2])
    ρ_AC = ptrace(ρ,[2,2,2],[2])
    ρ_BC = ptrace(ρ,[2,2,2],[1])
    
    S_ABC = von_neumann_entropy(ρ)
    S_AC = von_neumann_entropy(ρ_AC)
    S_BC = von_neumann_entropy(ρ_BC)
    S_C = von_neumann_entropy(ρ_C)
     I  = (S_AC + S_BC - S_C - S_ABC)
    return I
end

# Range of p values
p_values = range(0, stop=0.99, length=1000)
QCMI_values = [QCMI(p) for p ∈ p_values]
#Plot the function versus p
plot(p_values, QCMI_values, label = "I(A:B|C)")

Running the code, the following plot is obtained : Quantum Conditional Mutual Information

However, it is not correct because my supervisor told me I should obtain a different qualitative plot :

Real QCMI

How can I fix this problem ? Thanks to anyone can help me