diff --git a/types/types_xcpot_libxc.F90 b/types/types_xcpot_libxc.F90
index 5f3e98ca1211a519aade714b2ce0e6348d524d9b..9d4ccc9f3e74875e4d81cd1dfcc94d1ae8073f0b 100644
--- a/types/types_xcpot_libxc.F90
+++ b/types/types_xcpot_libxc.F90
@@ -111,9 +111,7 @@ CONTAINS
     ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0
     ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
     IF (xcpot%is_gga()) THEN
-       CALL judft_error("libxc GGA not implemented yet")
        IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
-       ALLOCATE(grad%vsigma,mold=grad%sigma)
        CALL xc_f03_gga_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,grad%vsigma)
        IF (xcpot%func_id_c>0) THEN
           ALLOCATE(vsigma,mold=grad%sigma)
@@ -170,10 +168,12 @@ CONTAINS
     INTEGER, INTENT (IN)         :: jspins,ngrid
     TYPE(t_gradients),INTENT(INOUT):: grad
     !For libxc we only need the sigma array...
-    IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr)
+    IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr,grad%laplace,grad%vsigma)
     ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
     ALLOCATE(grad%gr(3,ngrid,jspins))
-    
+    ALLOCATE(grad%laplace(ngrid,jspins))
+    ALLOCATE(grad%vsigma,mold=grad%sigma)
+
   END SUBROUTINE xcpot_alloc_gradients
 
 
diff --git a/xc-pot/libxc_postprocess_gga.f90 b/xc-pot/libxc_postprocess_gga.f90
index 30d604d01538496af5e9ca8c3adf8bfcedafec79..75a0001688648d8345d434185e16805d56687ca0 100644
--- a/xc-pot/libxc_postprocess_gga.f90
+++ b/xc-pot/libxc_postprocess_gga.f90
@@ -25,11 +25,11 @@ CONTAINS
     
     n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise
     nsp=SIZE(v_xc,1) !no of points
-    ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:atoms%lmax(n),n_sigma))
+    ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:sphhar%nlhd,n_sigma))
     vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
-    CALL mt_from_grid(atoms,sphhar,nsp,n,n_sigma,vsigma,vsigma_mt)
+    CALL mt_from_grid(atoms,sphhar,nsp/atoms%jmtd,n,n_sigma,vsigma,vsigma_mt)
     ALLOCATE(grad_sigma%gr(3,nsp,n_sigma))
-    CALL mt_to_grid(atoms,sphhar,vsigma_mt,nsp,n_sigma,n,.TRUE.,grad=grad_sigma)
+    CALL mt_to_grid(atoms,sphhar,vsigma_mt,nsp/atoms%jmtd,n_sigma,n,.TRUE.,grad=grad_sigma)
     
     CALL libxc_postprocess_gga(vsigma,grad,grad_sigma,v_xc)
   END SUBROUTINE libxc_postprocess_gga_mt
@@ -56,7 +56,7 @@ CONTAINS
     CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g)
     
     ALLOCATE(grad_sigma%gr(3,nsp,n_sigma))
-    CALL pw_to_grid(xcpot,3,.false.,stars,cell,vsigma_g,grad_sigma)
+    CALL pw_to_grid(xcpot,n_sigma,.false.,stars,cell,vsigma_g,grad_sigma)
     
     CALL libxc_postprocess_gga(vsigma,grad,grad_sigma,v_xc)
 
@@ -72,7 +72,7 @@ CONTAINS
     INTEGER:: i
     IF (SIZE(v_xc,2)==1) THEN !Single spin
        DO i=1,SIZE(v_xc,1) !loop over points
-          v_xc(i,1)=v_xc(i,1)-dot_PRODUCT(grad_sigma%gr(:,i,1),grad%gr(:,i,1))-vsigma(i,1)*grad%laplace(i,1)
+          v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_sigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1))
        ENDDO
     ELSE  !two spins
        DO i=1,SIZE(v_xc,1) !loop over points