Monthly Archives: September 2009

Natural Bond Orbitals (NBO) visualization

This is the second post on a series which will try to address common technical questions in computational chemistry that recursively appear on the CCL.

The Natural Bond Orbitals analysis is a powerful tool in population analysis calculations which is more robust than the traditional Mulliken approach, if for no better reason because its almost insensitive to the change of basis set while Mulliken’s P.A. is highly sensitive to basis set effects. Another advantage of the NBO analysis is that it provides a localized depiction of the electron density over a molecule, making it more related to chemists intuition. So far I have only worked with Gaussian 98, Gaussian 03 and only recently with Gaussian 09 in calculating NBO’s although it is possible also to perform them with GAMESS and the standalone NBO5.0 program created by Frank Weinhold. Visualizing them, however is never a straightforward process, and quite often we see more questions on the CCL than answers trying to address the matter. Most of the answers are concerned with what visualization programs to use but they seldom provide step by step instructions, furthermore most manuals are a bit cryptic about the procedure to plot this orbitals.

In Molekel 4.3

Make sure that the route section in your input file includes the following options: pop=nboread; gfoldprint (in case you are using G03 or even G09) or gfprint (in case you are still running G98), e.g.:

#P opt rhf/6-31G(d,p) pop=(full,nboread) gfoldprint geom=connectivity

This calculation requests a geometry optimization followed by a Natural Bond Orbital population analysis (with keywords to be read) using the Restricted Hartree-Fock method with the split valence basis set 6-31G(d,p)

NOTE: I have my own issues and concerns about the use of DFT along with NBO but maybe that will be treated in another post.

At the end of the input file, after the infamous blank line type the following:


You may include other keywords such as BNDIDX which generates a Wiberg bond index (order) matrix; or BOAO which generates the same matrix but in the Natural Atomic Orbitals basis. The PLOT option will generate a series of files with numerical extensions. The one you want to pay attention to is filename.47

If you are using Gaussian 03 or Gaussian 09 and STILL want to use Molekel 4.3 then you probably already know you have to change 03 for 98 on the header of the output file:

****************…              ****************
Gaussian 03:  x86…  –>       Gaussian 98:  x86-
2-Feb-2009                               2-Feb-2009
***************…                ****************

That aside, load your output file (filename.out or filename.log) on Molekel 4.3 as usual. Then go to Load -> nbo orb and load filename.47. Now, go to Compute -> Orbital and now select the Natural Bond Orbital you are interested in. This should do it! As an interesting exercise try computing the same orbitals (lets say HOMO and HOMO-1) with and without loading filename.47 in order to observe the difference between the shape of the MO’s and NBO’s. Molekel 4.3, though, is filled with bugs that will make it close unexpectedly, specially when running under Windows Vista. Sometimes the window closes because some sort of resolution problem, specially when taking snapshots (interestingly enough this happened to me when the background color was set to white) decrease the resolution of your monitor before taking the snapshot to prevent this problem. Some people complain about the look and feel of the latest molekel version so they stick to this old bugged one, so that is why I’m posting this method.

In Gasusview 3.0

This is the gaussview version I currently work with. When performing the NBO analysis on Gaussian (by the way, Gaussian 09 cannot be visualized with gaussview 3.x) use the savenbo option in the rout section, for example:

#P opt HF/6-31G(d) freq=noraman pop=(ful,nboread,savenbo) geom=connectivity

This will save the NBO coefficients into the checkpoint file. Load the output file with gaussview normally and then click on the Molecular Orbital icon (or go to Edit -> MO’s). This will open a new window with four tabs at the bottom of the molecule image. Click on the ‘New’ tab and load your checkpoint file. Gaussview will automatically format the chk file (this can cause some troubles when crossing architectures or platforms, so be careful to generate the proper formatted chk file!) Once load select the orbital you need to visualize and go to the ‘Visualize’ tab and click update; the orbital displayed is the Natural Bond Orbital. In this tab you can also adjust certain parameters like the isovalue (which is set to 0.02 by default) or the cube grid which controls how smooth the surface looks. UPDATE: According to John Keller from Alaska you can use this methodology when calculating with G09 and using Gview3.x, this also allows to visualize vibrations when this same software combination is used.

As usual, this post will be updated whenever I find some more useful information about the matter. Rate this post or leave a comment, just to know if you found it useful. Thanks!

Emergency talk. Calixarenes as guest molecules

My boss just told me a few days ago he may not be able to make it to a workshop in Bucharest, which is actually more of a bilateral conference between Romania and South Korea. So it seems I’m due up for making a presentation for this meeting on next Wednesday (that is the day after tomorrow!). The down side? I’m not available tomorrow (Sep. 15th) because I’m going to Bucharest to attend an Independence Day celebration thrown by the Mexican Embassy. Anyway, I had a layout of a presentation about the research on calixarenes I’ve done all year long., so I’m just filling in the slides with relevant data and some info about our facilities, as a way to also make some advertisement of our very own Babes-Bolyai University.

This talk will let me assess how much progress have we made so far and how much work do we still have to perform. Once again, quoting prof. Raymundo Cea-Olivares: Projects aren’t finished, they are dropped! It’s already the middle of September so it’s about time to wrap it up so we may have some published papers by the end of the year.  In the mean time, I better hurry up so I can celebrate Mexico’s independence 199th anniversary without any worries or due work.

Viva Mexico!

Polarizable Continuum Model (PCM) in G03

This is my first post on a series I have in mind regarding frequent questions on the CCL regarding the use of some computational chemistry software, mostly Gaussian. Readers are still encouraged to contact the Gaussian Help Desk for further (and more accurate) help.

Gaussian 03, the popular electronic structure calculation suite of programs, includes the necessary modules for performing calculations in a solvated environment using the continuum models approximations. Among such models, the Polarizable Continuum Model (PCM) is one of the most widely used methods since it meets a good compromise between accuracy and computation time. Nevertheless, Gaussian may not be the best option for performing such calculations (as opposed to other programs as COSMO) but it still can be very useful when used properly. Unfortunately there is a lack of specific info in the literature regarding the usage of the different variables involved in the cavity generation for G03; the newest version, G09, includes some improvements on the corresponding codes making PCM calculations more achievable. While browsing the CCL archives, it  is common to find more questions than answers and usually the same questions are posted over and over by different users over time. This post will get updated as needed.

I hope with this post I can summarize most of the common problems found in Gaussian regarding implicit solvation calculations as well as their respective solutions. Some of the solutions come from Gaussian technical support itself, so my best advise is always to address your questions directly to them. Keywords are typed in capital letters, variables in italics.

Brief background

Implicit solvent calculations imply the generation of a vacuum cavity inside a continuous and homogeneous dielectric field. The simplest model to do this is Onsager’s in which the molecule is treated as a dipole inside a spherical cavity (SCRF=DIPOLE in Gaussian use along the VOLUME keyword to generate the optimum radius for such cavity.) PCM calculations generate a cavity that relates more closely to the molecule’s shape by placing spheres on each atom or groups of atoms. Check the Gaussian link at the bottom of this post for further info; this is a troubleshooting post, not a tutorial on PCM.

Some common errors and their solutions

In order to get a better definition of a cavity it has been recommended to use the option SCRF=(READ,model,SOLVENT=solvent) with the following parameters to be read at the end of the input file:



Additionally we may include a third line indicating the kind of radii to be used on each atom to generate a sphere around it, the default option is Radii=UA0 (Topological United Atoms model) which treats functional groups as a single sphere. Including this line with Radii=UFF; Radii=Pauling or Radii=Bondi will treat each atom independently, which is very useful to use when some H atoms lye outside the UAKS sphere. The error message associated with this problem is: “Error message, treat H atom explicitly” see below

-> BldSpC: Error generating genealogic tree for sphere 309 at level 15

According to Gaussian’s Help Desk, this is a numerical error in the generation of the cavity. The use of fewer spheres (implicit H atoms for instance) is recomended, so if you are using RADII=PAULING or BONDI, delete that line. It is also recommended to use the NOSYMM keyword on the route section. This problem seems to have been addressed in G09.

-> Too many tesserae.  Increase the MxTs.

Try using the TSNUM keyword in the route section as SCRF=(TSNUM=num,…) This will modify the number of tiles to describe each sphere that makes the cavity.

-> AdVTs1: ISph=  500 is engulfed by JSph=  501 but Ae(  500) is not yet zero! Error in link301

Generation of cavity fails. Try using a different radii model (RADII=…) and/or the NOSYMMCAV keyword at the end of the file, via the SCRF=(Read,…) option. Also using the NOSYMM keyword in the route section can work. Once again using the OFAC=0.8 and RMIN=0.5 parameters is useful.

->  UA0: Hydrogen   40 is unbound. Keep it explicit at all point on the …

-> UA0: potential energy surface to get meaningful results.

The location of a certain H atom (number 40 in this case) lies outside the cavity placed on a functional group, so it must be treated explicitly by either changing the RADII= model or by placing a sphere on that particular atom alone through the SPHEREONH=40 (40 for this example) option via SCRF=(Read,…)

Additional remarks and suggestions

  • The use of spheres on functional groups is suggested for calculating energies, but for geometry optimization the use of a more sophisticated model in generating the cavity is encouraged.
  • Always pay attention to the value of the density lying outside the cavity, i.e. inside the dielectric. In G03 this value is labeled as “error on total polarization charges”. As a rule of thumb this value should be less than 0.05 for the calculation to be acceptable.
  • Just in case you are using a very old version of Gaussian, be aware that the keyword COSMORS doesn’t launch a COSMO-RS calculation (thermodynamics of solutes and solvents) but a CPCM calculation in a format that can be post-processed by COSMO software.
  • PCM calculations are highly parametrized so it’s useful to always have an experimental reference to which you can validate your choices in each calculation.

If you found interesting or helpful information in this post, please leave a comment however short. This will encourage me to keep gathering and posting this kind of information which in turn may be of help for other users, thanks.


%d bloggers like this: