problem with total charge density and graphite calculations

[Andrzej]

Dear users,

I have recently encountered some problem with calculations of my
crystal structure consisted of vdW binded layers, so I performed some
test for graphite, as at least I know what should be expected in this
case. There were single point calculations with a set of values of "c"
axis constant ruling the separation between graphite planes.

Reference results obtained by Sauer: J Comp. Chem, 29, (2008), 2088;
PBE+D; with a k-point mesh
equilibrium c-lattice constant: 6.4A, binding energy/layer/atom: 1.91
kcal/mol

my results:
simple cell (111)
- equilibrium c-lattice constant: 8.0 A, binding energy: -0.35 kcal/
mol (underestimated interaction)
- no minimum if pure PBE is used without D correction. (as expected)

cell doubled in all directions (2,2,2)
- equilibrium c-lattice constant (in reference to 111 cell) : 5.5 A,
binding energy: -4.14 kcal/mol (strongly overrestimated interaction)
- minimum on a curve generated using pure PBE at 6A  (in reference to
111 cell) with binding energy ca -1.2kcal/mol.

Can anybody explain me, why there is so large difference between 111
and 222 results? And why there is a minimum on the curve generated
using pure PBE. There should be no attractive interaction between
layers at this level, and curve should be totally repulsive.
Another  question, here are Mulliken charges from a run for 222 cell
with c constant 12A (6A in relation to basic cell):
     1  C     0.268132
     2  C    -0.404667
     3  C    -0.386350
     4  C     0.117564
     5  C    -0.032852
     6  C     0.549709
     7  C     0.143119
     8  C    -0.254631
     9  C     0.257086
    10  C    -0.403768
    11  C    -0.406458
    12  C     0.142254
    13  C    -0.020586
    14  C     0.537671
    15  C     0.137096
    16  C    -0.243319
    17  C     0.257094
    18  C    -0.403765
    19  C    -0.406460
    20  C     0.137081
    21  C    -0.020588
    22  C     0.537679
    23  C     0.142255
    24  C    -0.243318
    25  C     0.258142
    26  C    -0.404884
    27  C    -0.407575
    28  C     0.143326
    29  C    -0.021620
    30  C     0.538853
    31  C     0.138166
    32  C    -0.244386
  Total       0.000000

It looks like the atoms which should be equivalent are not seen as
equivalent during wave function optimization. I think that it may be
related to the mentioned above false interaction between carbon layers
at pure PBE level.

And last question,
when I performed similar calculations for 333 cell I got more or less
smooth curve, but for some separations total energy is dramatically
decreased and I observe a problem with total charge density counting.
Here are some results.:
c-axis constant for supercell(and in relation to basic cell); total
energy; total charge density;

15 (5); -613.4487600; 71e-10
18 (6); -614.160089 ; -1e-10
21 (7); -613.989957; -10e-10
24 (8); -623.06921; -0.00045 !!!
27 (9); -613.71077; +4.6e-6
45 (15); -613.31918 +11e-10

Why at some separations there is a sudden problem with charge density
counting?

I have played with GAPW, cutoff, restarting WFN obtained for other
separations, multigrid options... but to no effect. There are always
some points with wrong total charge density.

Here is one of my input for 333 cell:

&GLOBAL
  PROJECT 80
  PRINT_LEVEL MEDIUM
  RUN_TYPE ENERGY_FORCE
  FLUSH_SHOULD_FLUSH
&END GLOBAL

&MOTION
  &GEO_OPT
    MAX_ITER 500
#    OPTIMIZER CG
 &END GEO_OPT
&END MOTION

&FORCE_EVAL

  METHOD Quickstep
  &DFT
    WFN_RESTART_FILE_NAME 80-RESTART.wfn
    &MGRID
      CUTOFF 300
    &END MGRID
    &QS
       EPS_DEFAULT 1.0E-12
       MAP_CONSISTENT
    &END QS
    &SCF
      MAX_SCF 500
      EPS_SCF 1.0E-6
      EPS_EIGVAL 1.0E-5
      SCF_GUESS RESTART
      &OT
        MINIMIZER CG
#        N_DIIS 7
     PRECONDITIONER  FULL_ALL
      &END OT
        &OUTER_SCF
        MAX_SCF 20
        &END OUTER_SCF
   &PRINT
          &RESTART_HISTORY OFF
          &END
          &RESTART
             &EACH
                QS_SCF 0
                MD 10
             &END
             ADD_LAST NUMERIC
          &END
    &END
    &END SCF
    &XC
   &XC
       &XC_GRID
          XC_DERIV SPLINE2
       &END XC_GRID
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
      &vdW_POTENTIAL
         DISPERSION_FUNCTIONAL PAIR_POTENTIAL
         &PAIR_POTENTIAL
            TYPE GRIMME
         &END PAIR_POTENTIAL
      &END vdW_POTENTIAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
    A 7.3919999599  0.0000000000  0.0000000000
    B -3.6959999800  6.4016597501  0.0000000000
    C 0.0000000000  0.0000000000  24.0
      PERIODIC XYZ
    &END CELL
 &COORD
SCALED
.....
&END COORD
    &KIND C
      BASIS_SET DZVP-GTH-PBE
      POTENTIAL GTH-PBE-q4
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

the whole inputs with atoms positions as well as outputs for 333 runs
I have attached as gra333.tar file

I would be very grateful if somebody advise me how to force a proper
density counting. I have the same problems with other structures which
I calculate, and which have larger cell.  I did some tests with cutoff
up to 400. Increasing cutoffs to very large values is just not
possible in my cace. Using GAPW leads to even worse results. Can
anybody tell me what else I can change in my input to improve electron
counting?

Andrzej