Thursday, April 21, 2016

VASP CHGCAR difference and 2D plane cut through it

This post content is a copy from my awk workshop "Case studies" section.

From time to time, I need to calculate electron charge difference that is tabulated on a regular 3D grid. Common examples in my field will be Gaussian .cube files or VASP CHGCAR files. One can find some scripts, tools and programs that can do this in one way or another. In my case, again, I need something slightly different and... Well, doing it with awk is so simple that I never use other tools but the one I will mention here. With small changes I am able to subtract 6 files at the same time, and since the files tend to be large, I keep them compressed.


Some advantages of the script are:

  • extremely small memory footprint (the calculation is done line by line)
  • it does not need to know anything about the format, it only expects that all files have the same number of grid points (in each direction)
  • the number of atoms in the CHGCAR could be different
#!/usr/bin/awk -f

BEGIN {
  cmd1= "bzcat  "ARGV[1];
  cmd2= "bzcat  "ARGV[2];
  cmd3= "bzcat  "ARGV[3];

  NF=1;  while(NF>0){ cmd2 |& getline;       } cmd2 |& getline;
  NF=1;  while(NF>0){ cmd3 |& getline;       } cmd3 |& getline;
  NF=1;  while(NF>0){ cmd1 |& getline; print } cmd1 |& getline; print;

  do {
    cmd1 |& getline; split($0,d1)
    cmd2 |& getline; split($0,d2)
    cmd3 |& getline; split($0,d3)
    for(i=1;i<=NF;i++) printf "%18.11e ", d1[i]-d2[i]-d3[i]
    print""
  } while ($1*1==$1)

}

cmd1, cmd2 and cmd3 are the commands that will read the compressed files.

NF=1; while(NF>0){ cmd2 |& getline; } cmd2 |& getline; will fast forward to the grid data. The lines in the first file will be used as first lines in the new file. The rest is just reading line by line and calculating the difference column by column until it reaches the end of the file.

Bellow is an example of a 2D charge difference plot, produced by Python script that I wrote for this purpose. The script is specifically written for the case of water molecules, but could be easily tuned for more generic purposes.
FIg. 6 in International Journal of Quantum Chemistry116, ( 2016), 67-80; DOI: 10.1002/qua.25022

No comments:

Post a Comment