using horton do hirshfeld charge analysis

Date Tags python

In order to do hirshfeld charge analysis, I generate charge density cube file from cp2k (with GPW method, GAPW method has a problem generate charge density cube file). I post here the steps I use horton and this cube file to do hirshfeld charge analysis (follow this link ).

Here is how I use it on killdevil.

module load anaconda
source activate my_root

First, I need to set up a database for proatoms. I am using cp2k for atom calculation. I need H, O, Na, Cl in my system. So I generate the input files like this.

ppot.inc ( XXX_000_00 XXX stands for element number, GTH-PBE-qX which is the name of the pseudopotential)

001_000_00 GTH-PBE-q0
008_000_00 GTH-PBE-q6
011_000_00 GTH-PBE-q9
017_000_00 GTH-PBE-q7

valence.inc ( NNN_PPP_MM NNN: element number. PPP: atomic population(#electon). MM: spin multiplicity(2S+1) 1s2 2p1: orbital occupations)

001_001_02 1s1
001_002_01 1s2

008_005_02 1s2 2p1
008_006_03 1s2 2p2
008_007_04 1s2 2p3
008_008_03 1s2 2p4
008_009_02 1s2 2p5
008_010_01 1s2 2p6

011_005_02 1s2 2p1
011_006_03 1s2 2p2
011_007_04 1s2 2p3
011_008_03 1s2 2p4
011_009_02 1s2 2p5
011_010_01 1s2 2p6
011_011_02 1s2 2p6 2s1
011_012_01 1s2 2p6 2s2
011_013_02 1s2 2p6 2s2 3p1

017_013_02 1s2 2p1
017_014_03 1s2 2p2
017_015_04 1s2 2p3
017_016_03 1s2 2p4
017_017_02 1s2 2p5
017_018_01 1s2 2p6

template.inp

&GLOBAL
  PROJECT ATOM
  PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
  ATOMIC_NUMBER ${number}
  ELECTRON_CONFIGURATION (${mult}) CORE ${line:valence.inc}
  CORE none

  MAX_ANGULAR_MOMENTUM 1
  &METHOD
     METHOD_TYPE UKS
     &XC
       &XC_FUNCTIONAL PBE
         &PBE
           PARAMETRIZATION REVPBE
         &END PBE
       &END XC_FUNCTIONAL
     &END XC
  &END METHOD
  &POTENTIAL
      PSEUDO_TYPE GTH
      POTENTIAL_FILE_NAME ../../PBE_PSEUDOPOTENTIALS
      POTENTIAL_NAME ${line:ppot.inc}
  &END POTENTIAL
  &PP_BASIS
      BASIS_SET_FILE_NAME ../../BASIS_MOLOPT
      BASIS_TYPE CONTRACTED_GTO
      BASIS_SET DZVP-MOLOPT-SR-GTH
  &END PP_BASIS
  &PRINT
    &BASIS_SET ON
    &END
    &ORBITALS ON
    &END
    &POTENTIAL ON
    &END
  &END
&END ATOM

After preparing these files. I do

horton-atomdb.py input cp2k 1,8,11,17 template.inp

This will generate several folders and a run_cp2k.sh file. In order to run it on hopper (which supercomputer I use for cp2k) I add several lines in the head of this file.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#!/bin/bash
#PBS -N myjob 
###PBS -q regular
#PBS -q debug
#PBS -l mppwidth=24
###PBS -l walltime=96:00:00
#PBS -l walltime=00:30:00
#PBS -j oe
#PBS -V

# Note: if you want to use an mpi-parallel CP2K binary, uncomment the following
# line and fill in the right binary and mpirun script:
cd $PBS_O_WORKDIR
module load cp2k/2.5.1
#CP2K_BIN="mpirun -n4 cp2k.popt"
CP2K_BIN="aprun -n 24 cp2k.popt"

Then I do

qsub run_cp2k.sh

After the calculations are done, I copy the folder back to killdevil and using horton again. This will generate a atom.h5 file for hirshfeld calculation.

horton-atomdb.py convert

Now, I copy atom.h5 file and the cube file from cp2k calculation to another folder. I also modify the cube file based on the 3.4.1 in the manual. The command I use for the real calculation is 

horton-cpart.py cube output.h5:group {h,hi,he} atoms.h5

h, hi, he stands for different partitioning schemes.

To make the result readable, there is a tool to convert h5 file to csv file.

horton-hdf2csv.py output.h5 output.csv

The information I got based on hirshfeld analysis are in the file of output.csv.