plotting the charge density difference from cube file

[Xun Wang]

Hello,

I am a new user to cp2k and I would like to know how I can plot the atomic density difference from the cube file.

In my calculation of He dimer, I added the following lines to my input file to get the cube file. 

 &E_DENSITY_CUBE

 FILENAME ./rut_density

 &END E_DENSITY_CUBE 

I would like to get a plot as figure 2 in the following literature: http://scitation.aip.org/content/aip/journal/jcp/117/24/10.1063/1.1522715

I appreciate your help me with this!

Xun

[S Ling]

This is not a CP2K related problem. 

According to my understanding, you should first perform calculations for each single He atom placed at positions in the dimer, and you print out electron density cube files of each He atom, and then you sum them up, let's say you have cube file 1. Then you perform a calculation for the He dimer, and you print out the electron density cube file 2 (you named it rut_density). Finally, you take the difference of the two cube files by subtracting one from the other.

To manipulate the cube files, e.g. summation or subtraction, you can use the cubecruncher tool supplied with CP2K, see $CP2K/cp2k/tools/cubecruncher (where $CP2K is the main directory of your CP2K source code). Alternatively, you can use the cubman utility of Gaussian 09 to perform the cube file manipulations.

SL

[Xun Wang]

Thanks! I have one follow-up question about this: How can I get the one-dimensional density difference plot from the cube file ?

[S Ling]

You can plot it using the volume slice representation within VMD, or you can extract the values (let's say along the He-He bond axis and in the XZ plane) from the cube file and plot the data using other softwares. The cube file contains electron density difference at each grid point, and you should be able to find the number of grid points along each X/Y/Z axis in the begining of your cube file. Then you just need to correlate the electron density difference at each grid point with the actual cartesian coordinates of that grid point, and with these information, you can plot the electron density difference in whichever way you prefer. Be careful of the units in the cube file.

SL

[Samuel Andermatt]

Once you have the density difference in a separate cubefile you can also use the cubecruncher with -1d_profile.

[Xun Wang]

Thanks! I was trying to use the option -1d_profile, after I got the density difference file He-delta.cube,  however the following error occurs and I have no idea how to solve this. It seems an error of reading the data file. I could not locate the error. Could you help me with this! Thanks!

$ cubecruncher.x -i He-delta.cube -o out.cube -1d_profile z 

 Reading from cube          He-delta.cube

 Writing to   cube          out.cube

forrtl: severe (59): list-directed I/O syntax error, unit -5, file Internal List-Directed Read

Image              PC                Routine            Line        Source             

cubecruncher.x     000000000044A2A7  Unknown               Unknown  Unknown

cubecruncher.x     000000000046148F  Unknown               Unknown  Unknown

cubecruncher.x     000000000045FF1E  Unknown               Unknown  Unknown

cubecruncher.x     000000000041DDF7  Unknown               Unknown  Unknown

cubecruncher.x     00000000004147C3  Unknown               Unknown  Unknown

cubecruncher.x     000000000040DF1E  Unknown               Unknown  Unknown

libc.so.6          00007F71CE535D5D  Unknown               Unknown  Unknown

cubecruncher.x     000000000040DE29  Unknown               Unknown  Unknown

Thanks!!