Large residual charges in VDD method
Andrew-S-Rosen opened this issue · 4 comments
Andrew-S-Rosen commented
Apologies for the drastic edit as I continue to diagnose the issue and refine my question.
For the attached files below, I am getting a sum of VDD charges that is -50. This is the number of pseudopotential electrons in the system, whereas I expected to get a sum that is ~0. Do you know why I am getting the total number of electrons returned?
CHGCAR.gz
AECCAR0.gz
AECCAR1.gz
AECCAR2.gz
Input:
crystal AECCAR1
load AECCAR1
load AECCAR2
load as "$1-$2"
reference 3
voronoi
Output:
_ _ _ ___
(_) | | (_) |__ \
___ _ __ _ | |_ _ ___ ) |
/ __| | '__| | | | __| | | / __| / /
| (__ | | | | | |_ | | | (__ / /_
\___| |_| |_| \__| |_| \___| |____|
* CRITIC2: analysis of real-space scalar fields in solids and molecules.
(c) 1996-2019 A. Otero-de-la-Roza, A. Martin-Pendas, V. Lua~na
Distributed under GNU GPL v.3 (see COPYING for details)
Bugs, requests, and rants: aoterodelaroza@gmail.com
Website: https://aoterodelaroza.github.io/critic2/
If you find this software useful, please cite:
A. Otero-de-la-Roza et al., Comput. Phys. Commun. 185 (2014) 1007-1018.
A. Otero-de-la-Roza et al., Comput. Phys. Commun. 180 (2009) 157-166.
+ critic2 (development), version 1.0 (git:1.1dev-215-g84692b2)
host: Linux-5.3.18-150300.59.43_11.0.51-cray_shasta_c
date: Sat 18 Jun 2022 16:55:49
compiled dat: /global/homes/r/rosen/software/critic2dir/share/critic2
datadir: /global/homes/r/rosen/software/critic2dir/share/critic2
dic file: /global/homes/r/rosen/software/critic2dir/share/critic2/cif/cif_core.dic
...was found?: T
spglib: 1.13.0
libxc: <unavailable>
CRITIC2--2022/6/18, 20:56:50.793
critic2:1> crystal AECCAR1
%% crystal AECCAR1
* Crystal structure
From: AECCAR1
Lattice parameters (bohr): 10.905781 10.905781 10.905781
Lattice parameters (ang): 5.771091 5.771091 5.771091
Lattice angles (degrees): 134.065 134.065 66.986
Empirical formula: In(3) Ag(1)
Number of non-equivalent atoms in the unit cell: 3
Number of atoms in the unit cell: 4
Number of atomic species: 2
Number of electrons (with zero atomic charge): 194
Molar mass (amu, per unit cell): 452.322
Cell volume (bohr^3): 658.80883
Cell volume (ang^3): 97.62539
Density (g/cm^3): 7.69368
Space group (H-M): I4/mmm (139)
+ List of atomic species:
# spc = atomic species. Z = atomic number. name = atomic name (symbol).
# Q = charge. ZPSP = pseudopotential charge.
# spc Z name Q ZPSP
1 49 In 0.0000 --
2 47 Ag 0.0000 --
+ List of non-equivalent atoms in the unit cell (cryst. coords.):
# at = complete list atomic ID. xyz = Cartesian coordinates. spc = atomic species.
# wyck = wyckoff position. name = atomic name (symbol). mult = multiplicity.
# Z = atomic number.
# nat x y z spc wyck name mult Z
1 0.7500000000 0.2500000000 0.5000000000 1 2d In 2 49
2 0.5000000000 0.5000000000 0.0000000000 1 1b In 1 49
3 0.0000000000 0.0000000000 0.0000000000 2 1a Ag 1 47
+ List of atoms in the unit cell (cryst. coords.):
# at = complete list atomic ID. xyz = Cartesian coordinates. spc = atomic species.
# name = atomic name (symbol). Z = atomic number. nat = non-equivalent atom id.
# mol = molecular fragment.
# at x y z spc name Z nat mol
1 0.7500000000 0.2500000000 0.5000000000 1 In 49 1 1
2 0.2500000000 0.7500000000 0.5000000000 1 In 49 1 1
3 0.5000000000 0.5000000000 0.0000000000 1 In 49 2 1
4 0.0000000000 0.0000000000 0.0000000000 2 Ag 47 3 1
+ Lattice vectors (bohr)
a: -4.2554969554 4.2554969554 9.0949192858
b: 4.2554969554 -4.2554969554 9.0949192858
c: 4.2554969554 4.2554969554 -9.0949192858
+ List of atoms in Cartesian coordinates (bohr):
# at = complete list atomic ID. xyz = Cartesian coordinates. spc = atomic species.
# name = atomic name (symbol). Z = atomic number. dnn = nearest-neighbor distance.
# nat = non-equivalent atom id. mol = molecular fragment.
# at x y z spc name Z dnn nat mol
1 0.0000000000 4.2554969554 4.5474596429 1 In 49 6.0182 1 1
2 4.2554969554 0.0000000000 4.5474596429 1 In 49 6.0182 1 1
3 0.0000000000 0.0000000000 9.0949192858 1 In 49 6.0182 2 1
4 0.0000000000 0.0000000000 0.0000000000 2 Ag 47 6.0182 3 1
+ List of symmetry operations (16):
Operation 1:
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Operation 2:
-1.000000 0.000000 0.000000 0.000000
0.000000 -1.000000 0.000000 0.000000
0.000000 0.000000 -1.000000 0.000000
Operation 3:
0.000000 1.000000 0.000000 0.000000
0.000000 1.000000 -1.000000 0.000000
-1.000000 1.000000 0.000000 0.000000
Operation 4:
0.000000 -1.000000 0.000000 0.000000
0.000000 -1.000000 1.000000 0.000000
1.000000 -1.000000 0.000000 0.000000
Operation 5:
0.000000 1.000000 -1.000000 0.000000
1.000000 0.000000 -1.000000 0.000000
0.000000 0.000000 -1.000000 0.000000
Operation 6:
0.000000 -1.000000 1.000000 0.000000
-1.000000 0.000000 1.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Operation 7:
1.000000 0.000000 -1.000000 0.000000
1.000000 0.000000 0.000000 0.000000
1.000000 -1.000000 0.000000 0.000000
Operation 8:
-1.000000 0.000000 1.000000 0.000000
-1.000000 0.000000 0.000000 0.000000
-1.000000 1.000000 0.000000 0.000000
Operation 9:
-1.000000 0.000000 0.000000 0.000000
-1.000000 0.000000 1.000000 0.000000
-1.000000 1.000000 0.000000 0.000000
Operation 10:
1.000000 0.000000 0.000000 0.000000
1.000000 0.000000 -1.000000 0.000000
1.000000 -1.000000 0.000000 0.000000
Operation 11:
0.000000 -1.000000 0.000000 0.000000
-1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 -1.000000 0.000000
Operation 12:
0.000000 1.000000 0.000000 0.000000
1.000000 0.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Operation 13:
0.000000 -1.000000 1.000000 0.000000
0.000000 -1.000000 0.000000 0.000000
1.000000 -1.000000 0.000000 0.000000
Operation 14:
0.000000 1.000000 -1.000000 0.000000
0.000000 1.000000 0.000000 0.000000
-1.000000 1.000000 0.000000 0.000000
Operation 15:
-1.000000 0.000000 1.000000 0.000000
0.000000 -1.000000 1.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Operation 16:
1.000000 0.000000 -1.000000 0.000000
0.000000 1.000000 -1.000000 0.000000
0.000000 0.000000 -1.000000 0.000000
+ List of symmetry operations in crystallographic notation:
# number: operation ... ## order of rotation [axis]; [axis in Cartesian]
1: x,y,z ## 1
2: -x,-y,-z ## -1
3: y,y-z,-x+y ## 4 [1,1,0]; [-0.000,0.000,1.000]
4: -y,-y+z,x-y ## -4 [-1,-1,0]; [0.000,-0.000,-1.000]
5: y-z,x-z,-z ## 2 [1,1,0]; [0.000,0.000,1.000]
6: -y+z,-x+z,z ## -1 [1,1,0]; [0.000,0.000,1.000]
7: x-z,x,x-y ## 4 [-1,-1,0]; [-0.000,-0.000,-1.000]
8: -x+z,-x,-x+y ## -4 [-1,-1,0]; [-0.000,-0.000,-1.000]
9: -x,-x+z,-x+y ## 2 [0,1,1]; [1.000,0.000,0.000]
10: x,x-z,x-y ## -1 [0,1,1]; [1.000,0.000,0.000]
11: -y,-x,-z ## 2 [1,-1,0]; [-0.707,0.707,0.000]
12: y,x,z ## -1 [-1,1,0]; [0.707,-0.707,0.000]
13: -y+z,-y,x-y ## 2 [1,0,1]; [0.000,1.000,0.000]
14: y-z,y,-x+y ## -1 [1,0,1]; [0.000,1.000,0.000]
15: -x+z,-y+z,z ## 2 [1,1,2]; [0.707,0.707,0.000]
16: x-z,y-z,-z ## -1 [1,1,2]; [0.707,0.707,0.000]
+ List of centering vectors (1):
Vector 1: 0.000000 0.000000 0.000000
+ Crystal symmetry information
Space group (Hermann-Mauguin): I4/mmm (number 139)
Point group (Hermann-Mauguin): 4/mmm
Point group (Schoenflies): D4h
Holohedry: tetragonal
Laue class: 4/mmm
+ Cartesian/crystallographic coordinate transformation matrices:
A = car to crys (xcrys = A * xcar, bohr^-1)
-0.0000000000 0.1174950905 0.0549757490
0.1174950905 0.0000000000 0.0549757490
0.1174950905 0.1174950905 0.0000000000
B = crys to car (xcar = B * xcrys, bohr)
-4.2554969554 4.2554969554 4.2554969554
4.2554969554 -4.2554969554 4.2554969554
9.0949192858 9.0949192858 -9.0949192858
G = metric tensor (B'*B, bohr^2)
118.9360654890 46.4990481407 -82.7175568148
46.4990481407 118.9360654890 -82.7175568148
-82.7175568148 -82.7175568148 118.9360654890
+ List of fragments in the system (1)
# Id = fragment ID. nat = number of atoms in fragment. C-o-m = center of mass (bohr).
# Discrete = is this fragment finite?
# Id nat Center of mass Discrete
1 4 0.380762 0.873079 0.000000 No
+ This is a 3D periodic structure.
+ Atomic environment
Number of atoms (reduced cell/environment): 4 / 2376
Radius of (unit cell/environment) circumscribed sphere (bohr): 10.1080 / 75.2869
Maximum interaction distance (bohr): 35.3912
Covering regions:
Total number of regions: 800 (10 10 8)
Minimum region ID: -5 -5 -4
Maximum region ID: 4 4 3
Region side (bohr): 10.9604
Transformation origin (bohr): 2.1277,-2.1277,4.5475
Search offsets: 3375
Maximum search offset: 7
Average number of atoms per region: 2.9700
Maximum number of atoms in a region: 14
+ Vertex of the WS cell in cryst. coords. (18)
# id = vertex ID. xyz = vertex cryst. coords. d = vertex distance to origin (bohr).
id x y z d (bohr)
1 -0.640536 0.359464 -0.000000 6.53859929
2 0.359464 -0.640536 -0.000000 6.53859929
3 -0.359464 0.640536 -0.000000 6.53859929
4 -0.750000 -0.250000 -0.500000 6.22805295
5 -0.640536 -0.640536 -1.000000 6.53859929
6 -0.250000 -0.750000 -0.500000 6.22805295
7 -0.250000 0.250000 0.500000 6.22805295
8 0.359464 0.359464 1.000000 6.53859929
9 0.250000 -0.250000 0.500000 6.22805295
10 -0.359464 -0.359464 -0.000000 6.53859929
11 0.640536 -0.359464 -0.000000 6.53859929
12 0.250000 0.750000 0.500000 6.22805295
13 0.640536 0.640536 1.000000 6.53859929
14 0.750000 0.250000 0.500000 6.22805295
15 -0.250000 0.250000 -0.500000 6.22805295
16 -0.359464 -0.359464 -1.000000 6.53859929
17 0.250000 -0.250000 -0.500000 6.22805295
18 0.359464 0.359464 -0.000000 6.53859929
+ Faces of the WS cell (12)
# Face ID: vertexID1 vertexID2 ...
1: 1 3 12 13 8 7
2: 1 3 15 16 5 4
3: 1 4 10 7
4: 2 6 10 9
5: 2 11 14 13 8 9
6: 2 11 17 16 5 6
7: 3 15 18 12
8: 4 10 6 5
9: 7 10 9 8
10: 11 14 18 17
11: 14 13 12 18
12: 15 16 17 18
+ Lattice vectors for the Wigner-Seitz neighbors
# FaceID: Voronoi lattice vector (cryst. coords.)
1: 0 1 1
2: -1 0 -1
3: -1 0 0
4: 0 -1 0
5: 1 0 1
6: 0 -1 -1
7: 0 1 0
8: -1 -1 -1
9: 0 0 1
10: 1 0 0
11: 1 1 1
12: 0 0 -1
+ Lattice vectors for the Delaunay reduced cell (cryst. coords.)
a: -1 0 -1
b: 0 1 1
c: 1 0 0
Delaunay reduced cell lengths: 8.510994 8.510994 10.905781
Delaunay reduced cell angles: 112.967 112.967 90.000
+ Is the cell orthogonal? F
+ Is the reduced cell orthogonal? F
* List of scalar fields
+ Field number 0
Name: <promolecular>
Source: <generated>
Type: promolecular
Atoms in the environment: 2376
Use core densities? F
Numerical derivatives? F
Nuclear CP signature: -3
Number of non-equivalent critical points: 3
Number of critical points in the unit cell: 4
Alias for this field (2): $0, $rho0
This is the REFERENCE field.
* List of integrable properties (3)
# Id Type Field Name
1 v 0 Volume
2 fval 0 Pop
3 lval 0 Lap
* List of additional properties at critical points (0)
* List of core and pseudopotential charges for each field
# id type core? ZPSP
0 promol no
critic2:2> load AECCAR1
%% load AECCAR1
+ Field number 1
Source: AECCAR1
Type: grid
Grid dimensions : 96 96 96
Max. length Voronoi-relevant vector (bohr): 0.11360
First elements... 1.299366491349E+00 2.596550312426E+00 2.116855177832E+00
Last elements... 2.535805610145E+00 3.453233234926E+00 2.845799283748E+00
Sum of elements... 9.926157310707E+00
Sum of squares of elements... 8.684739290520E+03
Cell integral (grid SUM) = 0.00739140
Min: -1.56138005E+00
Average: 1.12193437E-05
Max: 6.73810221E+00
Interpolation mode (1=nearest,2=linear,3=spline,4=tricubic,5=smoothrho): 4
Use core densities? F
Numerical derivatives? F
Nuclear CP signature: -3
Number of non-equivalent critical points: 3
Number of critical points in the unit cell: 4
Alias for this field (1): $1
* Field number 1 is now REFERENCE.
critic2:3> load AECCAR2
%% load AECCAR2
load as "$+ Field number 2
Source: AECCAR2
Type: grid
Grid dimensions : 96 96 96
Max. length Voronoi-relevant vector (bohr): 0.11360
First elements... 2.730216181208E+00 3.771336493417E+00 3.216419679882E+00
Last elements... 3.629537692306E+00 4.547317924082E+00 3.996790863819E+00
Sum of elements... 6.715642382010E+04
Sum of squares of elements... 3.781930217824E+04
Cell integral (grid SUM) = 50.00728437
Min: -3.39070057E-02
Average: 7.59056078E-02
Max: 7.51720152E+00
Interpolation mode (1=nearest,2=linear,3=spline,4=tricubic,5=smoothrho): 4
Use core densities? F
Numerical derivatives? F
Nuclear CP signature: -3
Number of non-equivalent critical points: 3
Number of critical points in the unit cell: 4
Alias for this field (1): $2
critic2:4> load as "$1-$2"
%% load as "$1-$2"
reference + Field number 3
Source: <generated>
Type: grid
Grid dimensions : 96 96 96
Max. length Voronoi-relevant vector (bohr): 0.11360
First elements... -1.430849689859E+00 -1.174786180991E+00 -1.099564502049E+00
Last elements... -1.093732082161E+00 -1.094084689155E+00 -1.150991580070E+00
Sum of elements... -6.714649766279E+04
Sum of squares of elements... 3.159199095731E+04
Cell integral (grid SUM) = -49.99989296
Min: -1.77212666E+00
Average: -7.58943885E-02
Max: -8.20415622E-03
Interpolation mode (1=nearest,2=linear,3=spline,4=tricubic,5=smoothrho): 4
Use core densities? F
Numerical derivatives? F
Nuclear CP signature: -3
Number of non-equivalent critical points: 3
Number of critical points in the unit cell: 4
Alias for this field (1): $3
critic2:5> reference 3
%% reference 3
* Field number 3 is now REFERENCE.
* List of integrable properties (3)
# Id Type Field Name
1 v 0 Volume
2 fval 3 Pop
3 lval 3 Lap
* Summary of critical points for the reference field
Non-equivalent critical points = 3
Cell critical points = 4
Topological class (n|b|r|c): 3(4) 0(0) 0(0) 0(0)
Morse sum: 4
critic2:6> voronoi
%% voronoi
* Voronoi integration
+ Integrating atomic properties
+ Integrated property (number 1): Volume
+ Integrated property (number 2): Pop
+ Integrated property (number 3): Lap
* List of integrated atomic properties
+ The "Label" entries will be used to identify the integrable property
in the tables of integrated atomic properties. The entries with
an "x" will not be integrated.
# id Label fid Field Additional
1 Volume -- Volume
2 Pop 3 Field (v)
3 Lap 3 Laplcn (v)
* Key for the interpretation of table headings
# Id = attractor identifier
# cp/at = critical point/atom from the complete CP/atom list
# (lvec) = lattice translation from CP/atom to the main cell's representative
# ncp/nat = critical point/atom from the non-equivalent CP/atom list
# Name = atomic name
# Z = atomic number
# Position = attractor position
# Mol = molecule identifier (see corresponding table above)
# Volume = volume of the reference field basins.
# Pop (population) = the reference field integrated in its own basins.
# Lap = Laplacian of the reference field integrated in the reference field basins.
# $xxx = field xxx integrated in the basins of the reference field.
* List of attractors integrated
# Id cp ncp Name Z mult Position (cryst.)
1 1 1 In 49 -- 0.7500000 0.2500000 0.5000000
2 2 1 In 49 -- 0.2500000 0.7500000 0.5000000
3 3 2 In 49 -- 0.5000000 0.5000000 0.0000000
4 4 3 Ag 47 -- 0.0000000 0.0000000 0.0000000
* Integrated atomic properties
# (See key above for interpretation of column headings.)
# Integrable properties 1 to 3
# Id cp ncp Name Z mult Volume Pop Lap
1 1 1 In 49 -- 1.65341107E+02 -1.29013654E+01 8.48815172E-02
2 2 1 In 49 -- 1.64134792E+02 -1.28828182E+01 1.03197073E-01
3 3 2 In 49 -- 1.65243559E+02 -1.28746735E+01 1.65702455E-02
4 4 3 Ag 47 -- 1.64089369E+02 -1.13410359E+01 -2.04648836E-01
--------------------------------------------------------------------------------
Sum 6.58808826E+02 -4.99998930E+01 -4.94881913E-14
The sum of YT charges computed using the following input also yields a number >> 50, which I found odd.
crystal CHGCAR
load AECCAR0
load AECCAR2
load as "$1+$2"
reference 3
load CHGCAR
integrable 4
yt
# Id cp ncp Name Z mult Volume Pop Lap $4
1 1 1 In 49 -- 1.66818425E+02 6.19349153E+01 -9.86630407E+01 1.29251853E+01
2 2 1 In 49 -- 1.63212787E+02 6.18726747E+01 1.40387225E+01 1.28653767E+01
3 3 2 In 49 -- 1.65513501E+02 6.18885616E+01 1.24495299E+01 1.28797681E+01
4 4 3 Ag 47 -- 1.63264114E+02 5.57027366E+01 7.21747883E+01 1.13296700E+01
------------------------------------------------------------------------------------------------
Sum 6.58808826E+02 2.41398888E+02 -3.36754624E-10 5.00000001E+01
Just as a sanity check, I tried the Henkelman group's Bader code and got what appears to be reasonable values in return.
aoterodelaroza commented
The sum of YT charges is 50, which I presume is the number of valence electrons in your system. You need to look at the column $4 for the result of integrating field 4 (the CHGCAR). The "Pop" column is the integral of the all-electron density which, as mentioned in the other thread, has a very large grid integration error (and therefore is meaningless). Assuming you went with Ag pseudos with ZVAL=13 and In with ZVAL=11, your charges, calculated as ZVAL minus the $4 column, are about 0.1 for Ag and about -0.3 for In, which is not wholly unreasonable.
As for VORONOI, loading the AECCAR1 gives:
Cell integral (grid SUM) = 0.00739140
so it integrates to zero over the unit cell. This is inconsistent with the VASP manual, and it may be a quirk (or a bug) in VASP, because it should integrate to a number that is (very roughly) close to 50. It is supposed to be a valence density. Can you try with different pseudos? Maybe these datasets are missing info for this kind of calculation. I also suggest trying a simpler system with Ag and In separately.
In the meantime, you may want to try using critic2's atomic densities:
crystal AECCAR0
zpsp Ag 11 In 13
load AECCAR2 id ae2
load as core sizeof ae2 id core_pro
load as promolecular sizeof ae2 id ae_pro
load as "$ae_pro - $core_pro" id val_pro
load as "$val_pro - $ae2" id defrho
reference defrho
voronoi
This is not ideal, though, because they were generated in an entirely different way and therefore are not exactly comparable to the AECCAR2. Note the 0.15 residual charge when integrating the defrho over the unit cell.
aoterodelaroza commented
Also, can you please share the VASP version you are using, and maybe a private e-mail with the INCAR, KPOINTS, and POTCAR?
Andrew-S-Rosen commented
Ah, thanks! You're correct about the YT charges. Regarding the Voronoi charges, it does appear that there is some sort of bug with the AECCAR1 file in a certain set of older VASP versions that is causing the integrated density to be ~number of electrons instead of ~0. I've confirmed it's not restricted to just the example here. It looks like the "VDD" charges might just be regular Voronoi charges (i.e. the density integrated over the Voronoi cell, not the deformation density) in this case.
I will send you the inputs in a follow-up email later this week. I'll also dig into the version numbers. I just have to do a bit of file-digging to get the additional info we need to troubleshoot since I only have the volumetric data available at this very moment.
Andrew-S-Rosen commented
Thanks for your help, @aoterodelaroza. With the clues you provided, I was able to pinpoint the issue. It is not a VASP version or pseudpotential issue. It's an input issue. Using ADDGRID=.TRUE. causes the aforementioned issue in the AECCAR1 file. When this keyword is removed, the AECCAR1 is appropriate. As such, I will close this issue.