def _virial(self, T):
        """Virial coefficient

        Parameters
        ----------
        T : float
            Temperature [K]

        Returns
        -------
        prop : dict
            Dictionary with residual adimensional helmholtz energy:
                * B: ∂fir/∂δ|δ->0
                * C: ∂²fir/∂δ²|δ->0
        """
        Tc = self._constants.get("Tref", self.Tc)
        tau = Tc/T
        B = C = 0
        delta = 1e-200

        # Polinomial terms
        nr1 = self._constants.get("nr1", [])
        d1 = self._constants.get("d1", [])
        t1 = self._constants.get("t1", [])
        for n, d, t in zip(nr1, d1, t1):
            B += n*d*delta**(d-1)*tau**t
            C += n*d*(d-1)*delta**(d-2)*tau**t

        # Exponential terms
        nr2 = self._constants.get("nr2", [])
        d2 = self._constants.get("d2", [])
        g2 = self._constants.get("gamma2", [])
        t2 = self._constants.get("t2", [])
        c2 = self._constants.get("c2", [])
        for n, d, g, t, c in zip(nr2, d2, g2, t2, c2):
            B += n*exp(-g*delta**c)*delta**(d-1)*tau**t*(d-g*c*delta**c)
            C += n*exp(-g*delta**c)*(delta**(d-2)*tau**t*(
                (d-g*c*delta**c)*(d-1-g*c*delta**c)-g**2*c**2*delta**c))

        # Gaussian terms
        nr3 = self._constants.get("nr3", [])
        d3 = self._constants.get("d3", [])
        t3 = self._constants.get("t3", [])
        a3 = self._constants.get("alfa3", [])
        e3 = self._constants.get("epsilon3", [])
        b3 = self._constants.get("beta3", [])
        g3 = self._constants.get("gamma3", [])
        for n, d, t, a, e, b, g in zip(nr3, d3, t3, a3, e3, b3, g3):
            B += n*delta**d*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    d/delta-2*a*(delta-e))
            C += n*tau**t*exp(-a*(delta-e)**2-b*(tau-g)**2)*(
                    -2*a*delta**d+4*a**2*delta**d*(
                        delta-e)**2-4*d*a*delta**2*(
                            delta-e)+d*2*delta)

        # Non analitic terms
        nr4 = self._constants.get("nr4", [])
        a4 = self._constants.get("a4", [])
        b4 = self._constants.get("b4", [])
        Ai = self._constants.get("A", [])
        Bi = self._constants.get("B", [])
        Ci = self._constants.get("C", [])
        Di = self._constants.get("D", [])
        bt4 = self._constants.get("beta4", [])
        for n, a, b, A, B, C, D, bt in zip(nr4, a4, b4, Ai, Bi, Ci, Di, bt4):
            Tita = (1-tau)+A*((delta-1)**2)**(0.5/bt)
            Delta = Tita**2+B*((delta-1)**2)**a
            Deltad = (delta-1)*(A*Tita*2/bt*((delta-1)**2)**(
                0.5/bt-1)+2*B*a*((delta-1)**2)**(a-1))
            Deltadd = Deltad/(delta-1) + (delta-1)**2*(
                4*B*a*(a-1)*((delta-1)**2)**(a-2) +
                2*A**2/bt**2*(((delta-1)**2)**(0.5/bt-1))**2 +
                A*Tita*4/bt*(0.5/bt-1)*((delta-1)**2)**(0.5/bt-2))
            DeltaBd = b*Delta**(b-1)*Deltad
            DeltaBdd = b*(Delta**(b-1)*Deltadd+(b-1)*Delta**(b-2)*Deltad**2)
            F = exp(-C*(delta-1)**2-D*(tau-1)**2)
            Fd = -2*C*F*(delta-1)
            Fdd = 2*C*F*(2*C*(delta-1)**2-1)

            B += n*(Delta**b*(F+delta*Fd)+DeltaBd*delta*F)
            C += n*(Delta**b*(2*Fd+delta*Fdd)+2*DeltaBd*(F+delta*Fd) +
                    DeltaBdd*delta*F)

        prop = {}
        prop["B"] = B
        prop["C"] = C
        return prop