Problem in deriving charges from ACF.dat

[Elisa]

Hi to everybody,

I am trying to derive the partial charges of part of my system.
This latter, the one I am treating, is a pocket formed by a K+ ion (i.e. atom #84) coordinated by 2 carbonyl groups (i.e. atoms #13, #75) and 1 hydroxyl group (#31).

1) First, I performed a single-point calculation using Gaussian 09 in order to derive the ESP charges of my pocket. Here the keywords I have adopted to run the calculation:
#P B3LYP/6-31G* Pop=(ESP,ReadRadii, FormCheck) iop(6/33=2)
(where Radii is 1.8 Ang).

2) Then, I generated the cubefile by using cubegen in this way:
~$ cubegen Density=SCF filename.fchk 100 h

3) Finally, I have used the script (Version 1.0 01/11/17) provided by the Henkelman's group to perform Bader Charge Analysis (BCA).
~$ ./bader filename.cube

However, I have noticed that the electron number drastically changes passing from Gaussian to BCA. Indeed, while in Gaussian I have observed 324 e- when I run Bader script they become ~228 e-.
Why this happens if Gaussian performed an all-electron calculation?

If I understood correctly, I should subtract the number of valence-shell electrons from the total number of electrons returned by BCA. (ACF.dat file is reported below)
In my case, for instance:
K+ (8 e- in valence shell): 8 - 6.047251= 1.952749, it should be the K+ charge, but it is evidently incorrect.
Similar results were obtained for the oxygen atoms coordinating the K+ (which are #13, #75 and #31)
What am I doing wrong?

I don't understand how to derive the partial atomic charges from the so-obtained ACF.dat file

I really hope you can help me.
Thanks in advace for any contribution.

Regards,

Elisa

ACF.dat:

# X Y Z CHARGE MIN DIST ATOMIC VOL
--------------------------------------------------------------------------------
1 5.773000 -0.068431 5.364755 5.324592 0.976371 109.851111
------------------------------------------------------------------------------------------
12 6.079010 -1.433818 0.918292 1.188503 0.133246 16.803599
13 5.411278 -0.986524 -1.264083 7.558726 1.417691 123.854110 [carbonyl]
------------------------------------------------------------------------------------------
30 1.835404 -3.753498 0.188601 0.904707 0.304300 34.428063
31 0.025623 -2.821100 -3.119585 5.628525 1.169210 119.653210 [hydroxyl]
------------------------------------------------------------------------------------------
74 0.235727 5.512684 -1.982431 1.276686 0.116835 16.320737
75 0.934216 3.365554 -2.545042 7.768356 1.428152 129.551882 [carbonyl]
------------------------------------------------------------------------------------------
83 8.134373 8.616764 0.084571 7.824935 1.435170 139.933415
84 3.051417 -0.010162 -5.038998 6.047251 1.811534 132.400768 [K+]
--------------------------------------------------------------------------------
VACUUM CHARGE: 0.9258
VACUUM VOLUME: 43838.3192
NUMBER OF ELECTRONS: 227.9922

[graeme]

I think that the issue here is that it is hard to accurately represent the total charge density on a regular grid. Integration errors at the cusps in the charge density at the atomic nuclei can be very significant and are likely the cause of the error in the total number of electrons.

What we do to solve this problem is to do the Bader partitioning on the total charge density but then integrate only the valence charge in the Bader volumes. If you can get Gaussian to print the valence density, then you can use a command like:

bader valence_charge.cube -ref total_charge.cube

You can also add the -vac off flag to make sure that all of the charge is assigned to the atoms and not considered part of the vacuum region.

A finer grid can also help, but you typically have to go to a ridiculously large file in order to accurately integrate the total charge. The valence charge is not so problematic.

[Elisa]

Thank you Professor for your prompt reply.
As you suggested, I got Gaussian to print out the valence density, generating a new cube file. The command used is:
cubegen 0 MO=Valence filename.fchk valence_charge.cube 200 h
Then I have used the command you proposed me:
bader valence_charge.cube -ref total_charge.cube -vac off
Unfortunately, I have encountered a problem and the following message:
GRID BASED BADER ANALYSIS (Version 1.0 01/11/17)

OPEN ...
-------------------------------------------
RUN TIME: 0.20 SECONDS
CALCULATING BADER CHARGE DISTRIBUTION
0 10 25 50 75 100
PERCENT DONE: **********************
REFINING AUTOMATICALLY
ITERATION: 1
EDGE POINTS: 222658
REASSIGNED POINTS: 23903
---------------------------------------------
ITERATION: 21
CHECKED POINTS: 35
REASSIGNED POINTS: 0
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
bader 00000000004D1DE5 Unknown Unknown Unknown
bader 00000000004CFA07 Unknown Unknown Unknown
bader 0000000000486094 Unknown Unknown Unknown
bader 0000000000485EA6 Unknown Unknown Unknown
bader 00000000004371B6 Unknown Unknown Unknown
bader 000000000043ADB0 Unknown Unknown Unknown
bader 00000000004F6280 Unknown Unknown Unknown
bader 000000000041175D Unknown Unknown Unknown
bader 0000000000400602 Unknown Unknown Unknown
bader 000000000040050E Unknown Unknown Unknown
bader 00000000004F7ABB Unknown Unknown Unknown
bader 00000000004003E9 Unknown Unknown Unknown

Reading the description of the options (‘bader -h’) I have noticed that the flag ‘-ref’ is used in order to process the VASP files.
Since in the website ‘CODE: BADER CHARGE ANALYSIS’ it is stated that both VASP CHGCAR files and Gaussian CUBE files can be processed by the program, how can I fix this problem?
Thanks for your help.

Regards,
Elisa

[graeme]

If you can post the charge density files (total and valence) I can try to find the reason for the error.

[Elisa]

I really appreciate your taking the time.
Here the charge density files.

Thank you really much,
Elisa

[graeme]

There appears to be a problem with your valence charge density. Notice that the file is much larger than the total charge. There are also extra lines below the atomic positions. I'm don't know enough to help with gaussian, but see if you can get a valence charge density in the same form as your total charge density. You can try loading it into a program such as vesta to make sure that it looks correct.

[Elisa]

Thank you so much for spending your time to help me.
While I was trying to solve the incoherence in charge densities, I drew two simpler systems (PDBs annexed as .txt), one system with K+ (sysA.txt - 12 atoms) and the other one without K+ (sysB.txt - 11 atoms). The system is composed by one water molecule and two formaldehyde molecules. I run a g09 singlepoint calculation for both systems:

> #P HF/6-31G Density=current Pop=ReadRadii (sysA)
> #P HF/6-31G Density=current (sysB)

Using formchk I generated the .fchk files necessary to cubegen. The command:

cubegen 0 Density=SCF sysA.fchk sysA.cube 100 h
cubegen 0 Density=SCF sysB.fchk sysB.cube 100 h

is started for both systems. Thus, the cube files obtained (annexed hereto) have been processed by bader program:

bader -i cube sysA.cube -vac off


>> GRID BASED BADER ANALYSIS (Version 1.0 01/11/17)
-------------------------------------------------------
CALCULATING BADER CHARGE DISTRIBUTION
-------------------------------------------------------
NUMBER OF BADER MAXIMA FOUND: 245024
SIGNIFICANT MAXIMA FOUND: 17
VACUUM CHARGE: 0.0000
NUMBER OF ELECTRONS: 31.95122

bader -i cube sysB.cube -vac off

>> GRID BASED BADER ANALYSIS (Version 1.0 01/11/17)
--------------------------------------------------------
CALCULATING BADER CHARGE DISTRIBUTION
--------------------------------------------------------
NUMBER OF BADER MAXIMA FOUND: 20
SIGNIFICANT MAXIMA FOUND: 20
VACUUM CHARGE: 0.0000
NUMBER OF ELECTRONS: 31.99057

The result obtained for 'sysB' is absolutely reliable (ACFsys.txt - section sysB), while the result obtained for 'sysA' (ACFsys.txt - section sysA) made me worry something was wrong with K+.
First of all, since the number of electrons in both systems is the same, I deduced that K+ has no valence electrons. Is it correct? Then, looking at the charge values I noticed that for 'sysB' they properly represent the chemical properties of the molecule, such as the C=O bond's polarity, while for 'sysA' the values are completely different (see ACFsys.txt).

Since I was not sure if the problem is due to bader analysis or Gaussian calculation, I have tried to calculate Merz-Singh-Kollman (MK) charges. The command:
> #P HF/6-31G Density=current Pop=(ReadRadii,MK) ; (sysA)

and the result is absolutely reliable and in line with bader charge analysis for 'sysB':

1 C 0.552165
2 O -0.603550
3 O -0.956013
4 H 0.473366
5 C 0.537151
6 O -0.601536
7 K 0.937164
8 H 0.047672
9 H 0.042916
10 H 0.489696
11 H 0.035983
12 H 0.044987

My question is whether it is a problem in generating the cube file, or in calculating charges with bader program.

Thank you really much.
Regards.

Elisa

[Elisa]

Hi,
I found a way to include all electrons (i.e., core and valence) of my system in the cubefile (using the flag “fdensity” on cubegen), not only the valence electrons. Indeed, so far the problem was the cubefile related to valence charge density and the incorrect charge of K+ obtained by running 'BCA code'; thus we tried to fix it by means of generating an “all-electrons” cubefile. We concluded that maybe the bad calculation was ascribable to the total number of electrons. Now I got that number (see upper).

Then I run the 'bader charge analysis code' obtaining the following ACF.dat file

-ACF.dat:

# X Y Z CHARGE MIN DIST ATOMIC VOL
-------------------------------------------------------------------------------- --------------
1 5.776898 -0.083283 5.363159 0.077185 0.000000 1.352014
2 5.259092 -1.918201 5.304578 0.574305 7.320362 39.208397
3 4.250011 0.969312 5.812914 0.464907 3.613017 29.599442
4 7.377497 0.071656 6.389280 1.345775 2.623376 98.890143
5 6.468546 0.693386 2.812029 0.392538 3.647716 7.436075
6 5.253470 2.242976 2.199758 0.285429 2.443969 9.464096
7 9.212439 1.587196 2.734550 0.298339 0.864000 5.070051
-----------------------------------------------------------------------------------------------
80 5.108032 8.486632 -3.665952 0.931310 0.458923 1444.787069
81 5.973518 7.770416 0.209876 4.677338 0.620169 87.151908
82 5.019202 7.411379 2.299913 9.006466 1.276941 802.315978
83 8.137264 8.611320 0.094603 8.879272 1.293874 1798.775181
84 3.051913 -0.007663 -5.041673 17.928845 2.050907 3706.652255 (K+ )
--------------------------------------------------------------------------------
VACUUM CHARGE: 0.0000
VACUUM VOLUME: 0.0000
NUMBER OF ELECTRONS: 323.3257

The charge on K+ is 17.928845 with the final net charge of 0.6743. This is clearly unphysical.

Moreover the total number of electrons should be 324 instead of 323.3257.
Hence, the problem still persists.
I am still in trouble but I would like to fix permanently the problem. Any suggestion?

Thanks a lot.
Regards,

Elisa

[graeme]

It is very hard to properly integrate the total charge density as represented on a grid. The cusps at the nuclei make this particularly hard. Our solution is to use the total charge density to do the Bader partitioning but then integrate the valence charge in those Bader volumes. The valence charge is much smoother and is easier to represent accurately on a grid. To do this, run a command such as:

bader charge_valence -ref charge_total