Working with cubes can be tedious. I need to show a change in electronic density of a MOF. For that I made two cubes for neutral and charged MOF. Then took their difference using cube_tools, like this.
import numpy as np
from cube_tools import cube
# Load the cube files using the cube class
cube1 = cube('mof_opt_0.0.cube')
cube2 = cube('mof_opt_2.0.cube')
# Subtract the data from cube1 from cube2
cube_diff = cube('mof_opt_2.0.cube')
cube_diff.data = cube2.data - cube1.data
# Get z-axis data and find indices where z > 13.3 (jellium density)
z_indices_above_threshold = np.where(cube_diff.Z > 13.3)[0]
# Remove densities above z = 13.3 by setting them to zero
for idx in z_indices_above_threshold:
cube_diff.data[:, :, idx] = 0
# Save the modified cube data to a new file
cube_diff.write_cube('cdd.cube')
Once I have got the charge density difference and opened it in VMD, I realised that one part of my MOF is right at the border of a periodic cell, so that part of density was split. So, I used a terminal command to shift the cube like this “cube_tools -t -24 -36 0 cdd.cube”. I had to shift manually positions of the atoms by taking into account the voxels size. Next challenge was hiding half of my MOF to clear the view. So I used this tcl syntax in VMD:
vmd > mol clipplane normal 0 0 0 {1 0 0}
vmd > mol clipplane center 0 0 0 {3 0 0}
vmd > mol clipplane status 0 0 0 1
vmd > mol clipplane normal 0 1 0 {1 0 0}
vmd > mol clipplane center 0 1 0 {3 0 0}
vmd > mol clipplane status 0 1 0 1
vmd > mol clipplane normal 0 2 0 {1 0 0}
vmd > mol clipplane center 0 2 0 {3 0 0}
vmd > mol clipplane status 0 2 0 1
Here is the result – density is almost homogeneously spread over my MOF upon charging.