DDEC6 Charges, Valence Density Mismatch

[Nicholas Laws]

Hi CP2K Community,

I am trying to obtain time-dependent DDEC6 charges for a mixed QM/MM simulation in which an EMI-BF4 ion pair (76 valence e, 26 core e) collides with a gold slab at various kinetic energies. The QM region is treated with ωB97X-V/aug-cc-pVTZ and all-electron potentials in GAPW. During production (input ESCollision.inp, attached) I print a cube with:

&PRINT

  &E_DENSITY_CUBE

    DENSITY_INCLUDE SOFT_DENSITY

    NGAUSS 0

    STRIDE 1 1 1

    FILENAME valence_density

  &END

&END

and a 5-level PW grid of:

&MGRID

  NGRIDS     5

  CUTOFF     800

  REL_CUTOFF 80

  COMMENSURATE .TRUE.

&END

The log shows:

Electronic density on regular grids:        -76.7658266487       25.2341733513

Core density on regular grids:              102.0000000000       -0.0000000000

Hard and soft densities (Lebedev):         -694.5370264067     -669.3028566617

Total Rho_soft + Rho1_hard - Rho1_soft     -101.9999963937

Total charge density (r-space):               0.0000036063

Total Rho_soft + Rho0_soft (g-space):         0.0000033449

Overlap energy of the core charge distribution:               0.00002215241243

Self energy of the core charge distribution:               -720.29674329623936

Core Hamiltonian energy:                                    -75.86850888340894

Hartree energy:                                             209.46847518535864

Exchange-correlation energy:                                -26.02788481604846

Hartree-Fock Exchange energy:                               -28.87249823090302

Dispersion energy:                                            0.25782644744755

so the cube integrates to ≈77 e instead of the expected 76 e. Chargemol therefore exits with

ncore =       26.0000

nvalence =       76.0000

pixelvolume =    1.3079E-03

numerically integrated valence density =    7.6767E+01

sum_valence_occupancy_correction =    0.0000E+00

checkme =    7.6659E-01

The electrons are not properly accounted for. Either the grid in your electron density input file is too coarse, you have specified the incorrect net charge in the chargemol_job.m file, or the number of core electrons for each atom has not been set up correctly.

Program will terminate.

My questions are:

1. Is there any keyword combination in GAPW that prints only the valence density (ρ_soft + ρ₁,hard − ρ₁,soft) so the integral matches 76 e, or must I run a GPW single-point for every MD frame?

2. Is there a published benchmark for the PW cutoff required to converge ωB97X-V + aug-cc-pVTZ?

Any hints, keywords, code pointers, or literature would be greatly appreciated!

Many thanks in advance,

Nick