def METHOD_NAME(self):
    n2c = mol.nao_2c()
    numpy.random.seed(10)
    ngrids = 100
    coords = numpy.random.random((ngrids,3))*2
    coords = coords[70:75]
    dm = numpy.random.random((n2c,n2c))
    dm = dm + dm.T.conj()
    aoLa, aoLb, aoSa, aoSb = r_numint.eval_ao(mol, coords, deriv=1)
    rho0a = numpy.einsum('pi,ij,pj->p', aoLa[0], dm, aoLa[0].conj())
    rho0b = numpy.einsum('pi,ij,pj->p', aoLb[0], dm, aoLb[0].conj())
    rho0 = rho0a + rho0b
    aoL = numpy.array([aoLa[0],aoLb[0]])
    m0 = numpy.einsum('api,ji,bpj,xab->xp', aoL.conj(), dm, aoL, lib.PauliMatrices)
    ni = r_numint.RNumInt()
    rho1 = ni.eval_rho(mol, (aoLa[0], aoLb[0]), dm, xctype='LDA')
    self.assertAlmostEqual(abs(rho1[0].imag).max(), 0, 9)
    self.assertAlmostEqual(abs(rho0-rho1[0]).max(), 0, 9)
    self.assertAlmostEqual(abs(m0 - rho1[1:4]).max(), 0, 9)