Question about methods requiring promolecular density

[Andrew-S-Rosen]

Hello @aoterodelaroza,

Thanks for making such a great code. I am a new user and interested in using Critic2 for some of the methods that involve a promolecular density. Specifically, I am most interested in calculating (Voronoi) deformation density charges and think that Critic2 might be the best way for me to get them. Deformation density is mentioned very briefly in passing in the manual, but it's not clear to me how to go about doing this. I am using a VASP output.

Would you be able to provide some additional details about how to calculate the deformation density for a periodic system?

[aoterodelaroza]

If you are not in a hurry I can probably write one or two examples about this in the next few days and stop being cryptic about it :wink: . It should be indeed possible to integrate the Voronoi cell using the structural variables, although if you want to use this very often it is probably a better idea if I just make a new keyword for it. The tools are all in the program already - wouldn't take too long.

[Andrew-S-Rosen]

Certainly, if you are willing to go through the effort to share some details about this, I am most happy to wait for however long that takes! Thank you!!

I would plan to be using this for quite a number of materials. A keyword would be fantastic! But even if that doesn't end up being feasible in the short-term, any additional insights would be great.

Looking forward to learning more!

[aoterodelaroza]

OK, I implemented the VORONOI and HIRSHFELD keywords. Also, an example of how to calculate VDDs specifically. Let me know if this works.

[Andrew-S-Rosen]

@aoterodelaroza, wow thank you so much! I am going to try it out right now. I'm so excited!

[Andrew-S-Rosen]

@aoterodelaroza -- I tried compiling critic2 and got the following error. This is with the intel 19.1.2.254 compiler. Any suggestions about how to resolve this one? Thanks in advance.

cd /global/homes/r/rosen/software/critic2/build/src && /opt/cray/pe/craype/2.7.10/bin/ftn -DHAVE_READLINE -I/global/homes/r/rosen/software/critic2/build/src/cubpack -I/global/homes/r/rosen/software/critic2/src/qhull -I/global/homes/r/rosen/software/critic2/build/src/json-fortran  -O3 -qopenmp   -c /global/homes/r/rosen/software/critic2/src/hirshfeld@proc.f90 -o CMakeFiles/critic2.dir/hirshfeld@proc.f90.o
/global/homes/r/rosen/software/critic2/src/hirshfeld@proc.f90(133): error #6645: The name of the module procedure conflicts with a name in the encompassing scoping unit.   [HIRSH_NOGRID]
  subroutine hirsh_nogrid()
-------------^
compilation aborted for /global/homes/r/rosen/software/critic2/src/hirshfeld@proc.f90 (code 1)

EDIT: It compiles perfectly fine with gfortran 11.1.2.0. You can disregard this comment as a result.

[aoterodelaroza]

Ooops. I pushed the fix for this. Please git pull and try again.

[Andrew-S-Rosen]

Thanks! I was able to get it compiled and ran the VDD charge analysis per the example! I'm going to go ahead and close the issue.

[Andrew-S-Rosen]

@aoterodelaroza -- actually, I have a few short clarification questions if you don't mind just to make sure I'm understanding it correctly.

1. Why do we use AECCAR1-AECCAR2 instead of AECCAR1-CHGCAR? Just by trying to draw parallels, in the Bader method it looks like we treat rho(r) as the CHGCAR, but here we are treating it as AECCAR2. Edit: I guess CHGCAR and AECCAR2 are pretty much the same in hindsight.
2. I understand why we are subtracting AECCAR2 from AECCAR1 in the integration since that's the (negative) deformation density. However, when calculating the Voronoi cell of atom A (i.e. where the integration is taking place), why do we use AECCAR1-AECCAR2 as the reference to do that instead of AECCAR0+AECCAR2 or CHGCAR?

[aoterodelaroza]

AECCAR2 and CHGCAR both contain the valence electron density. CHGCAR is the smooth pseudo-density actually used in the plane-wave calculation and AECCAR2 is the "true" valence electron density reconstructed with the PAW transformation. See this - AECCAR2 would be the equation at the bottom and CHGCAR only the first term on the rhs of that equation. The AECCAR2 has very sharp peaks at the nuclear positions (like the true all-electron density), the CHGCAR does not.

The thing with grids is that whatever field you integrate on them needs to be smooth - otherwise the quadrature error will be very high. So if you want to integrate the valence number of electrons, as in the Bader charges, you need to use the CHGCAR. For the VDD, you need to integrate the deformation density, and the promolecular valence density written by VASP in the AECCAR1 is comparable to the AECCAR2, i.e. it is not pseudized. Therefore, you need to subtract AECCAR1 - AECCAR2. The deformation density calculated in this way is (allegedly) smooth because the sharp peaks in AECCAR1 cancel with those in AECCAR2 when you take the difference, because the core density is mostly unaffected by the chemical environment of an atom. However, you should check the integral of the deformation density over the whole cell, just in case (what I called the "residual sum" up above).

[Andrew-S-Rosen]

Thanks for the super clear explanation. That all makes perfect sense to me!