Hirschfeld-I Analysis

The goal here is to implement the Hirschfeld-I (iterative) analysis described by Bultinck et al (paper1) (paper 2). The only application to solid state calculations was by Laennarts et al

Things needed to do this:

  1. We have to be able to break up each atom into a separate "type"; I have written a code wsx:/homes/jschrier/hirscheld-i/separate_atomtypes.f90 that can do this. Takes one argument, which is the density file name. Basically this parses an Abinit version >=5.7 density file and gives each atom a unique type, following the description of the density file format. This also generates a file containing the "basenames (i.e., atom names as characters)" for each of the types of atoms, that will be used to select the correct all-electron density file later on. This code is run with the separate_atomtypes.r file, found in the same directory, and using the following command: ./separate_atomtypes.r mz3X8o_DEN
  2. New strategy required for atomic density: Because LDA/GGA do not yield stable cations (but we need them!) we have resorted to a new strategy, using a slightly modified version of C. Froese Fischer's atomic Hartree Fock code (Computer Physics Communications, 43 (1987) 355) so that it outputs the formatted variables.
    1. Code is in /homes/jschrier/atomic/hf_js.r
    2. The easiest way to run this code is to feed in an input file. See HF input files
    3. Run ./hf_js.r <input*

From before:Josh's Notes

  1. Run perl /homes/jschrier/hirscheld-i/run_hi.pl to perform hirshfeld-i analysis
    1. Have all the ae files in the current directory
    2. Run perl /homes/jschrier/hirshfeld-i/interp_charges.pl mz318o_DEN.basenames -init to generate the initial (all neutral) ae files
    3. perl run_hi.pl mz318o_DEN.sep mz318o_DEN.basenames
  2. The above can be tedious if you want to do a lot of iterations, so we recommend using a "meta-script" that will perform the whole Hirshfeld-I process.
    1. perl /homes/jschrier/hirshfeld-i/interp_charges.pl mz318o_DEN.basenames -init (initialize system, as before)
    2. perl meta_run.pl mz318o_DEN.sep mz318o_DEN.basenames (will run 100 iterations) (note that you need to have copies of both meta_run.pl and run_hi.pl in the current working directory where you are analyzing your charge. (also please note that if you are missing some of the required ae files, then it will start spewing out lots of error messages. Kill the job using ctrl-c if that happens)
    3. Visualize output files, e.g., xmgrace -nxy charge_iterations.dat (there is also a file of charge_delta_iterations.dat which just shows the difference between the previous iteration. The x-axis is the iteration number, and each atom gets its own line in these plots)
    4. In case of charge sloshing, it is helpful to restart from the average of the LAST 50 iterations. You can do this with perl damp_oscillation.pl mz318o_DEN.basenames charge_iterations.dat (this generates both an output file of average_charges.dat, as well as the new ae files to be used to restart the Hirshfeld-I process)
    5. Then re-run the Hirshfeld-I process, perl meta_run.pl mz318o_DEN.sep mz318o_DEN.basenames
    6. When you are happy with the result, plot the mz318o_DEN.basenames charge_iterations.dat (e.g., using xmgrace -nxy mz318o_DEN.basenames charge_iterations.dat) to see if the charges have converged; usually they do, in which case you can just read them out from the tail of the follow
  3. In some very early calculations, there were oscillations of the charges, and so we extracted the average charges. perl final_avg.pl mz318o_DEN.basenames charge_iterations.dat will extract the average of the FIRST 20 charge iterations, and write it to final_average_charges.dat

If error messages stating that a certain .ae file does not exist, but it has been made, remember that UNIX is case-sensitive, so check for incorrect capitalization.

Hirshfeld-I evaluation using VASP

(this refers to codes that are in the carver:/project/projectdirs/nanoHC/hirshfeld-i directory)

  1. Perform your VASP calculation using PAW. To output the core and valence densities (on the NGXF, NGYF, NGZF grid), use the LAECHG = .TRUE. ** option in your INCAR file.
  2. There is a problem in that the core densities fluctuate quickly, so the default grid will not integrate to the correct total charge. To fix the integral problem, perform a calculation where you set **ISTART = 2, ICHARG = 0 (this will reload from your existing WAVECAR file and use that to generate the initial charge density; this saves you from having to do an entire SCF calculation again, and then set NGXF = 999, NGYF=999, NGZF=999 (don't really set these to 999 as you will certainly run out of memory. Instead, double the default NGXF, etc. values that are listed in the OUTCAR file. My experience with KVO3 is that the default setting gives ~0.8% too much charge, and doubling the default NGXF settings (when in PREC = medium mode) reduces this to ~0.3%. These files rapidly get larger, but if you can afford it, make these as large as possible. Note that in order to reload the wavefunctions, you must keep the planewave energy fixed and use the same number of processors. Also, it is a good idea to specify NBAND by hand in these restarts to keep it the same as in your initial calculation.
  3. This will output AECCAR0 and AECCAR2 files. Combine these using perl ../chgsum.pl AECCAR0 AECCAR2 (this code is from the Graeme Henkelman's group @UT Austin) This will produce a file CHGCAR_sum that contains the sum of the valence and core densities
  4. Make a new directory, and move your CHGCAR_sum file to it. Then change its name to CHGCAR (be careful about overwriting the existing CHGCAR file if you care about it). Then convert this file to abinit format using parchg2den.r (this code is from Steve Erwin@NRL). Note that this utility program only works with VASP 4.x outputs, and will die if used with VASP 5.2 outputs. Not a problem, as the default vasp that is loaded on nersc is 4.x. This will produce a file named DEN containing the density
  5. Use separate_atomtypes.r.paw to generate the sep-files. (For now this only works with KVO3, we haven't generalized to other atom types yet when dealing with VASP-generated inputs; this executable will also work on AbInit generated inputs as above.) Note that the density files generated by parchg2den.r do not contain atom information, so you have to create the basenames file by hand. This is the same format as the basenames file that get produced automatically frm the AbInit-generated inputs: first line is the number of atoms, then a blank (comment) line, then the names of the atoms on each successive line, in lowercase, on each successive line of the file. E.g., ../../../../hirshfeld-i/separate_atomtypes.r.paw DEN will generate a DEN.sep file and assumes there is a DEN.basenames file present (note that the code will give you a warning to remind you of this
  6. From here, you can simply proceed with the usual Hirshfeld-I calculation discussed above.

(VASP H-I charges do not currently work)

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License