from qutip import *
from pylab import *
%pylab inline

B0=(tensor(basis(2,0),basis(2,0))+tensor(basis(2,1),basis(2,1))).unit()

def negat(rho):
    xx=partial_transpose(rho,[0,1])
    vv=xx.eigenenergies()
    n=0.0
    for i1 in range(0,rho.shape[0]):
        n=n+abs(vv[i1])-vv[i1]
    return n

plist=linspace(0,1,50)
nlist=zeros(len(plist))

for i1 in range(len(plist)):
    p=plist[i1]
    rho=(1.0-p)*ket2dm(B0)+p/4.0*tensor(qeye(2),qeye(2))
    nlist[i1]=negat(rho)

plot(plist,nlist,'-',label='', linewidth=4)
xlabel('p',fontsize=20)
ylabel('N',fontsize=20)
show()