Printing Core Charges for Bader analysis

[Satish Kumar]

Hello CP2K users

I was trying to obtain "proper?" Bader charges (Using Henkelman's program) for a CO2 molecule. Right now, I am getting the charge associated with C atom to be zero. In VASP, we can add the LAECHG=TRUE tag that will print the core charges too (http://theory.cm.utexas.edu/henkelman/code/bader/). Doing so in VASP gives a reasonable (around 2 electrons) charge associated with C atom. After looking up the CP2K Google group I tried a couple of other options in an attempt to get the total (valence + core) charge density. I tried the following print sections:

1)

&PRINT

    &E_DENSITY_CUBE

      STRIDE 1

    &END E_DENSITY_CUBE

&END PRINT

2)

  &PRINT

   &E_DENSITY_CUBE

      STRIDE 1

      TOTAL_DENSITY .TRUE.

    &END E_DENSITY_CUBE

  &END PRINT

3) This one gave wrong results I think (number of electrons estimated from Bader was 0). 

  &PRINT

    &TOT_DENSITY_CUBE

      STRIDE 1

    &END TOT_DENSITY_CUBE

  &END PRINT

Input file:

&FORCE_EVAL

  METHOD Quickstep

  &DFT

    UKS

    BASIS_SET_FILE_NAME /home/siyemperumal/GTH_BASIS_SETS_5-12-10

    POTENTIAL_FILE_NAME /home/siyemperumal/GTH_POTENTIALS_5-12-10

    WFN_RESTART_FILE_NAME co2-RESTART.wfn

    &MGRID

      CUTOFF 1600

    &END MGRID

    &QS

      WF_INTERPOLATION ASPC

      EXTRAPOLATION_ORDER 3

    &END QS

    &SCF

     EPS_SCF 1.E-6

     SCF_GUESS RESTART

     MAX_SCF 500

     &OT T

       PRECONDITIONER FULL_SINGLE_INVERSE

       MINIMIZER DIIS

       LINESEARCH 3PNT

     &END OT

    &END SCF

    &XC

      &XC_FUNCTIONAL PBE

      &END XC_FUNCTIONAL

      &VDW_POTENTIAL

         DISPERSION_FUNCTIONAL PAIR_POTENTIAL

         &PAIR_POTENTIAL

            TYPE DFTD3(BJ)

            PARAMETER_FILE_NAME /home/siyemperumal/Research/cp2k/package/cp2k/data/dftd3.dat 

            REFERENCE_FUNCTIONAL PBE

    &PRINT_DFTD MEDIUM

    &END PRINT_DFTD

         &END PAIR_POTENTIAL

      &END VDW_POTENTIAL

    &END XC

  &PRINT

#Tried the above three combinations

  &END PRINT

  &END DFT

  &SUBSYS

    &CELL

      ABC 12.00 13.00 14.00 

    &END CELL

    &TOPOLOGY

    &END TOPOLOGY

    &COORD

C         0.0000153980        0.0000010425        0.0006891492

O         1.1762907217        0.0000004039        0.0006107401

O        -1.1763057650       -0.0000013603        0.0006092107

    &END COORD

    &KIND O

      BASIS_SET DZVP-MOLOPT-GTH 

      POTENTIAL GTH-PBE-q6

    &END KIND

    &KIND C

      BASIS_SET DZVP-MOLOPT-GTH 

      POTENTIAL GTH-PBE-q4

    &END KIND

  &END SUBSYS

&END FORCE_EVAL

&GLOBAL

  PROJECT co2

  RUN_TYPE ENERGY 

  #RUN_TYPE GEO_OPT 

  PRINT_LEVEL LOW

&END GLOBAL

 &MOTION

  &GEO_OPT

    MAX_ITER 200 

    MAX_FORCE 0.0009725  #0.05 eV/A

    OPTIMIZER BFGS 

  &END GEO_OPT

 &END MOTION

Related Output that is printed:

  Electronic density on regular grids:        -15.9999999999        0.0000000001

  Core density on regular grids:               15.9999999994       -0.0000000006

  Total charge density on r-space grids:       -0.0000000005

  Total charge density g-space grids:          -0.0000000005

Thank you for any ideas. 

[Matt W]

Hi,

I've never tried this.  I think you would have to run a GAPW calculation (all electron) to get core electrons explicitly into the calculation. The you would have to use your second type of print statement with total_density added.

Pseudopotential calculations don't have a model of the frozen core density like PAWs do, so you can't get this from a GPW calculation easily (You could maybe hack it somehow by changing the sign of the nuclear charge or something in a TOTAL_DENSITY_CUBE run).

Matt 

[Satish Kumar]

Thank you for the idea. Will dig deeper into this.

--
You received this message because you are subscribed to a topic in the Google Groups "cp2k" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cp2k/q75YaO3Z4ro/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cp2k+uns...@googlegroups.com.
To post to this group, send email to cp...@googlegroups.com.
Visit this group at https://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/d/optout.

[S Ling]

Hi

You may try this tool, see: 

https://sourceforge.net/projects/ddec/ 

Although the tool was originally developed to perform partial charge analysis using another method (see http://pubs.acs.org/doi/abs/10.1021/ct100125x), the tool can produce a cube file which contains the core electron density based on some reference densities (included in the DDEC program). You can then perform Bader charge analysis on this new cube file using the tool developed by Henkelman. Please have a look at the manual included in the DDEC program. You may need to change some of the parameters (number of core electrons) in the DDEC program, as the GTH pseudopotentials used by CP2K may correspond to different numbers of core electrons.

SL

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.

[Satish Kumar]

Will try it out. Thank you.

You received this message because you are subscribed to a topic in the Google Groups "cp2k" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cp2k/q75YaO3Z4ro/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cp2k+uns...@googlegroups.com.

[David T]

Hi
I downloaded the code but the documentation is a bit encrypted.
I can use the code for calculating the charges but I can't write any reconstructed cube file.

Do you know how to do it? Can you please explain it to me?

Cheers
Davide

[Satish Kumar]

I do not think there are any options in the job_control.txt file to write the reconstructed cube files in DDEC6 program. I do not know if you want to tweak the program to make it write these reconstructed files. If you have a CP2K valence electron density file with you, some relevant files may be module_format_valence_cube_density.f08 or module_add_missing_core_density.f08 (not sure)?

I believe there will a paper soon by Prof. Manz, which will provide you with the core densities that you can use to generate your own reconstructed cube files. 

On a different note for interested users, when CP2K files are generated with large cutoff values and if you end up getting Fortran runtime error/End of file error, the default variable value of read_buffer_size in module_global_parameters may need to be increased from its default value. 

HTH.

--
You received this message because you are subscribed to a topic in the Google Groups "cp2k" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cp2k/q75YaO3Z4ro/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cp2k+uns...@googlegroups.com.

[S Ling]

Hi Davide,

I looked at this a long time ago when it was still version 3. As suggested by Satish, it looks like you will need to hack the DDEC code in order to produce a new cube file with re-inserted core density. Maybe the best way is to get in touch with Thomas Manz who may already have a solution for that.

Alternatively, if you are working with porous materials (e.g. MOFs), you may try the REPEAT scheme by Campana et al (DOI: 10.1021/ct9003405) which we recently implemented into CP2K. To use the method, just insert the following section into the &FORCE_EVAL section of your input:

   &PROPERTIES

     &RESP

       USE_REPEAT_METHOD  T

       &SPHERE_SAMPLING

         AUTO_VDW_RADII_TABLE  UFF

         AUTO_RMIN_SCALE     1.0

       &END SPHERE_SAMPLING

       &PRINT

         &PROGRAM_RUN_INFO  SILENT

         &END PROGRAM_RUN_INFO

       &END PRINT

     &END RESP

   &END PROPERTIES

See https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP.html#USE_REPEAT_METHOD for more information.

SL