This short article will guide you through the creation of images like the one above based on charge density calculations using Gaussian09. For visualization, I will use Jmol. The principal steps are:
- Get a conformation of the molecule you are interested in.
- Optimize this geometry until you find a local minimum and calculate the charge density distribution at that point.
- Generate a formatted Gaussian checkfile.
- Generate a cubefile.
- Use Jmol for rendering.
Get a conformation
As you work with a specific molecule, I’m sure that you have plenty of conformations at your hands. In principle, there are numerous ways how to obtain the necessary Gaussian input file, but I like to use vmd’s force field toolkit for this purpose. In order to start this way, you have to have a PDB file. Fire up vmd and select “Extensions”->”Modeling”->”Force Field Toolkit”. Select the tab “Opt. Geometry” and define both input and output files. A click on “Write Gaussian Input File” will do the obvious. (Of course, writing the input file by hand is a possibility, as well. In particular, as PDB and GAU are quite similar formats. But I personally tend to forget about the linebreaks…)
For testing purposes, you may use this Gaussian input file (which is for hexyfluoroethane)
%chk=basecheck.chk %nproc=8 %mem=4GB # MP2/6-31G(d) Opt=(Redundant) SCF=Tight Geom=PrintInputOrient Density=MP2 Pop=Full <qmtool> simtype="Geometry optimization" </qmtool> 0 1 C1 -0.23199999332427979 0.7289999723434448 0.0 C2 0.23199999332427979 -0.7289999723434448 0.0 F3 0.23199999332427979 1.3519999980926514 1.0920000076293945 F4 0.23199999332427979 1.3519999980926514 -1.0920000076293945 F5 -1.5709999799728394 0.777999997138977 0.0 F6 -0.23199999332427979 -1.3519999980926514 -1.0920000076293945 F7 -0.23199999332427979 -1.3519999980926514 1.0920000076293945 F8 1.5709999799728394 -0.777999997138977 0.0
Optimize and calculate
Now Gaussian should work on this input file for about a minute. The “Opt=Redundant” keyword does optimizations in internal coordinates (see Gaussian manual for details). “SCF=Tight” sets a strict convergence criterion for the geometry optimization and “Geom=PrintInputOrient” is not that hard to understand. The more interesting keywords are the last two: “Density=MP2” is responsible for calculating the charge density from the MP2 simulation data. This is set just to be sure that you actually use a compatible representation (and not i.e. SCF). Of course, both “MP2/6-31G(d)” and “Density=MP2” can be adjusted according to the manual to fit your requirements. The last keyword “Pop=Full” tells Gaussian to dump all population related data to the checkpoint file which is defined in the first line of the sample input file. Otherwise Gaussian will not retain this information.
Generate a formatted checkfile
As the cubegen tool of the Gaussian suite works on formatted checkfiles, we have to convert the simulation run output. In general, it is a good idea to convert everything to formatted checkfiles as the normal ones are sensible to architecture or version changes.
formchk -3 basecheck.chk base.fchk
does the trick.
Generate a cubefile
Now that you have all necessary data, it has to be processed such that Jmol (or other tools) are able to read them. The formatted checkfile contains a lot of information (depends on the route description you used in the Gaussian input file). For this purpose, Gaussian comes with the cubegen tool which extracts information from the formatted check file:
cubegen 0 density=MP2 base.fchk basemp2.cube 0 h
Rendering
Now start Jmol and select “File”->”Console”. After
load basemp2.cube
you should see the molecule only. With
isosurface plane (atomno=1) (atomno=2) (atomno=3) color range 0 .2 "basemp2.cube"
you should get something similar to the image on top of this page. The three atom numbers define a plane (here, it’s the only plane of hexafluoroethane that contains four atoms) more conveniently than it could be done without the console of Jmol.