import pylab as pl
import numpy as np
from gpaw.setup_data import SetupData
from gpaw.xc_functional import XCFunctional
from gpaw.utilities import unpack, unpack2
from gpaw.grid_descriptor import AERadialGridDescriptor as AERGD


np.set_printoptions(precision=4, suppress=True)
xc = XCFunctional('PBE')

s0 = SetupData('H', 'PBE')

s = s0.build(xc, 2, 1, None)
Delta0_ii = unpack(s.Delta_pL[:, 0].copy())
print Delta0_ii

#for phit_g, phi_g in zip(s0.phit_jg, s0.phi_jg):
phit = s0.phit_jg[0]
phi = s0.phi_jg[0]

nt = phit**2
n = phi**2
rgd = AERGD(s0.beta, s0.ng)

r_g = rgd.r_g

pl.plot(r_g, n - nt, 'r')
#pl.plot(r_g, nt)
#pl.show()
#sdkfjsdkjf

print 'n - nt', (n - nt).sum()

print 'sum', (nt - n).sum()
print 'integral', np.dot(n - nt, rgd.r_g**2 * rgd.dr_g), len(rgd.r_g)
import pylab as pl
pl.plot(rgd.r_g, (n - nt) * rgd.r_g**(2) * rgd.dr_g)
#pl.plot(np.arange(len(rgd.r_g)), rgd.r_g, 'ro')
pl.show()

#r_g = s0.rgd.r_g
#dr_g = s0.rgd.dr_g
integral = (( (1./4./np.pi)**.5 * (phi**2 - phit**2)) * rgd.r_g**2 * rgd.dr_g).sum()
#print integral