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

http://molekel.cscs.ch/wiki/pmwiki.php

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:

$NBO PLOT $END

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 propper 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.

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:

OFac=0.8

RMin=0.5

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.

References

http://www.ccl.net

http://www.gaussian.com/g_tech/g_ur/k_scrf.htm

http://www.cosmologic.de

Back to Romania

After two months of working at Pécsi Tudomanyegyetem in the research group of Prof. Kunsagi Sandor, I go back now to my previous office at the Babes-Bolyai University in Cluj-Napoca, Romania. Working here at Pécs has been a stimulating experience as well as a very productive one. We managed to obtain enough data as to prepare a couple of papers, which should be ready for submission in about a month. Several calculations on the hosting properties of a selected series of calixarenes were performed in order to further shed light on their skeletal properties with and without guests. NBO calculations were also carried out to assess the bonding properties within these molecules. Most calculations were carried out at the HF/6-31G(d) level of theory and only some DFT (B3LYP) were performed for comparison only. We expect to further understand the mechanisms by which calixarenes work as molecular recognition agents.

Be this a public recognition of the Hungarian hospitality, specially that of Prof. Kunsagi Sandor.

Kösönöm!

Chemistry and comic books

Science permeates into the collective mind of a society not only through school but also through the form of popular media such as the TV or in the case of this post comic books.

For a long time now, the group of  John Selegue and James Holler at the University of Kentucky have a website named as the comic book periodic table of chemical elements. In this clickable periodic table we can browse scans of the pages of different comic books in which the corresponding element is mentioned.

Clicking on each element will display some options for different comic books related to it but for every comic book only the page in which the corresponding element is referenced, is shown, therefore one is not able to read the entire comic book. Nevertheless, for the hardcore comic book fan there are, in many cases, insights to what the page displays putting the comic book, as well as the chemistry within it, in context.

Through the use of the physical properties of chemical elements, many comic book writers have created characters that base their identities in such properties. Also in some occasions, the chemical knowledge of a character helps to the story development. In both cases, chemical concepts are being -literally- illustrated and, ultimately imprinted in the collective mind.

In some cases it seem that there was a deliberate attempt to create a character with a storyline that revolved around the properties of the corresponding chemical element. Such is the case of all the “Metal Men” series in which a group of individuals have super powers related to the main characteristics of each’s corresponding metal. In other cases, such as in mainstream comic books like Batman or Superman, the inclusion of chemical knowledge is brought in by the supervillian who is usually a “mad scientist” trying to take over the world. This vision of scientists with power to enslave the human race probably arose from the atomic era as a consequence of rapid weapon development having the Manhattan project as an imediate antecedent.

The use of popular art forms has the benefit of reaching a larger audience and hence it also has the responsability of not distorting scientific facts into pseudo scientific ones.

Once again this post comes from my memories from the chemistry faculty back at UNAM and the classes of Dr. Raymundo Cea-Olivares who introduced it to us.

Explaining Entropy can be a mess…

Another scientific concept that is hard to grasp by laypeople and that to my opinion has been the center of much distortion in the chemistry classroom, is the thermodynamical function Entropy, S.

More often than not, S is said to be a measure of “disorder” and people just take it! If one was to define disorder then one would have to also define order: Is my apartment too entropic? what about my life? Does nature understand order in the same way as we do? How do we understand order inside a living cell where many molecules and organelles are floating around? If indeed S was a measure of disorder then, why is it important to measure it?

Entropy in a nutshell. There have been many attempts to define S in a way young students may understand it, yet tracing parallelisms with ordinary every-day-life concepts is hard and often leads to miss conceptions.  A student of mine once asked: “if entropy is always increasing, how come bodies tend to cool down?” he meant to ask how come the translation motions of a molecular ensamble tended to decrease (and with this achieving “order”.)

Prof. Mayo Martínez-Kahn at UNAM in Mexico wrote a very interesting paper about Entropy in the local journal of the Chemistry School, “Educación Química”. The paper was entitled “The tombs of Entropy” as a reference to the widely known fact that in Boltzmann’s tomb his famous equation relating Entropy to the partition function Q, is engraved. Prof. Martínez then ventures in imagining how would other tombs from people who have made contributions to the concept and notion of S would look like. I remember distinctively the one of Sadi Carnot’s in which his famous thermodynamic cycle was displayed.

Entropy in so many words is a function that describes how many different energy levels are available in a thermodynamic system. The more levels, the higher the entropy. It also describes the spontaneity of a process to occur since in nature a system always tends to undergo changes that increase its entropy along with that of its surroundings.

How come Gibbs’ free energy or Helmholtz don’t cause such confusions? my guess is because nobody has attached an every-day-word to them!

PS. It is still important to make scientific concepts permeate into the general audience. Recently decesead comedian George Carlin mentioned Entropy in the following video…

Are theoretical/computational chemists real chemists?

Although classifications and rankings do not define the range or scope of work a person/group can do, they are usually popular to look at, unfortunately these views are sometimes misleading in the decision making process, e.g. funding. How expensive is in average our science compared to that of say organic chemists? It indeed depends on the computer facilities at hand and how much are they used, as well as the kind of synthesis performed by this hypothetical org. chem.

The line between theoretical and computational is becoming diffuse, and the branches are increasing (chemoinformatics, chemometrics, bioinformatics and so on.) But since they are somewhat intangible, the results of our calculations seem to often fall into the black box category and hence dismissed or taken with reserves.

There is no doubt that a lot of companies are making a lot of money by selling all kinds of software related to chemistry, from electronic structure calculations to molecular sorting and recognition as well as data analysis. It cannot be denied the enormous amount of scientific effort put in develop such products. Also, pharmaceutical companies make a full use of computational methods as a way to find new drugs in shorter times.

Not too long ago, in the Computational Chemistry List, there was a discussion about whether or not computational chemists were ‘arm chair scientists’ and I’ve been meaning to post on the subject since then. Are we arm chair scientist? maybe so. Is that a bad thing? Of course not! My father is a retired mathematician and I know for a fact that he seldom left his desk, let alone his office. The idea of being an ‘arm chair scientist’ was just obvious for him and there was no need to ever question it. A friend of mine from Mexico told me that there are very few experimental physicists there because for years it was cheaper to make theoretical physics than experimental. I think nowadays students are encouraged to undertake experimental physics in order to balance the ratio.

From my point of view computational chemists can be divided into the following three major categories:

1) Those who develop new theories regarding the electronic structure of molecules as well as new models for modeling its behavior and properties.

2) Those who develop new codes for solving long standing problems. Within this branch we can find also those who develop new basis sets, density functionals, pseudopotentials, semi-empirical methods, force fields, etc.

3) Those (such as myself) who use the products of the previous two branches and apply them into solving problems of chemical interest.

We certainly don’t tackle chemistry the same way our experimentalist colleagues but then again they don’t tackle chemistry the way they used to a hundred years ago. I remember that prof. Raymundo Cea-Olivares at UNAM used to say that chemistry has become ephemeral, since synthesis lasts a day and then the job turns into physics when acquiring and analyzing spectra in order to assess the compound’s identity.

One thing I’ve noticed is that apparently it has become trendy for every research group to develop its own code and those are not good news for people with good knowledge of theoretical chemistry as well as good chemical insight but with poor programming skills.

Just pouring my opinions on the subject…

Inorganic Chemistry paper

A new paper is coming out! These are always good news for someone whose productivity is evaluated by the number of his/her publications, and in my case the pleasure is double. It turns out that despite the fact a scientist is continuously working, there isn’t always the possibility of having results put “out there”. Back in Mexico we have the “National Researchers System” in which the National Council for Science and Technology encourages us to keep on working by providing economical stimuli and evaluating our productivity by, yes, the number of papers published each year. For three years I worked for a private research center in which fundamental science was tackled as far as it was economically possible and cost effective. Practically all of the work carried out there was confidential since the company it belonged to is leader in its market, not only in Mexico but in Latin America and southern USA!  At this facility some papers were published from time to time (most of them came from research in molecular modeling) but not without struggle against the administrators. We can thank for most of the struggle (and the papers!) to Dr. Armando Gama-Goicochea, a great physicist as well as a great friend of mine.

Anyway, the bottom line here is that I’m excited about having a new paper coming out again, even if I’m nowhere close to being first author. This three year paper fastening seems to be over and let us hope it’s only the first of various now that I’m holding a posdoc position here in Romania.

In this paper we tackled the bonding properties of some Aluminum complexes with three chalcogeno triazoles. The electrostatic potential mapped onto a density surface of one of those compounds is currently shown in the header of this blog, btw. We concluded that the bonding in such compounds is mainly covalent as opposed to the more conventional electrostatic notion prevailing for such hard atoms. In order to get this information we resorted to Natural Bonding Analysis calculations with the RHF method and somewhat large basis sets in order to get a full description of the electronic density.

I very much like these systems in which several bonding possibilities occur. The fact that nature is chosing one out of many has always a reason which can be assessed by our models and may serve us to learn how to modify it’s behaviour.

 

Coordination Diversity of Aluminum Centers Molded by Triazole Based Chalcogen Ligands

Inorg. Chem., 2009, 48 (13), pp 5874–5883

Jocelyn Alcántara-García, Vojtech Jancik, Joaquín Barroso, Sandra Hidalgo-Bonilla, Raymundo Cea-Olivares, Rubén A. Toscano and Mónica Moya-Cabrera

Moving out

As part of the research project I’m working on here at the Babes-Bolyai University in Romania, I must now pay a visit to prof. Sandor Kunsagi-Mate in Pecs, Hungary. I’ve met prof. Kunsagi in several occasions before and I’m looking forward to working with him in his research team, although unfortunately he will have to leave next week to the far east to check on some of his multiple collaborations around the world. I’m also looking forward to use their computer facilities which I think are physically located in Budapest. Living near Budapest is something I’m also very much excited about since it’s one of my favorite cities in the world.
I haven’t read much about the University I’m going to, other than it’s the oldest university in Hungary dating back from the days in which King Matei Corvin ruled the Hungarian Empire. King Mathias was born in Cluj-Napoca, by the way. That’s a good 600 years now!
It’s always exciting to start anew in a completely different environment with a new research team. This time my chances of learning a bit of the language are pretty slim. In fact I’ve heard a lot of Hungarian here in Cluj and I haven’t learned too much, other than some isolated words.
Despite all the problems we face on a daily basis here at UBB I will miss Romania while I’m gone; it could be for a month or two. Romania has crawled under my skin, so even if I -often- complain about their ways and customs, it is now a part of me and I think I’m starting to feel homesick.
Pecs will be a nice adventure, I’m sure.

Wheel? I think knot!

Once again an awful title. This post follows my previous one on graphs and chemistry, and it addresses an old idea which I have shared in the past with many patient people willing to listen to my ramblings.

It is a common conception/place to state that the wheel was the invention that made mankind spring from its more hominid ancestors into the incipient species that would eventually become homo sapiens; that it was the wheel, like no other prehistoric invention or discovery, what made mankind to rise from its primitive stage. I’ve always believed that even if the wheel was fundamental in the development of mankind, man first had to build tools to make wheels out of something; otherwise they would have been just a good theoretical conception.

But even despite the fact that building tools was in itself a pretty damn good start, I strongly believe that mankind’s first groundbreaking invention were knots. For even a wheel was a bit useless until it was tied to something. From my perspective, the invention of the wheel was an event bound to happen since there are many round shaped things in nature: from the sun and the moon to some fruits and our own eyes. Achieving the mental maturity of taking a string (or a resembling equivalent of those days) and tie it, whether around itself or to something, was, in my opinion, the moment in which the opposable thumbs of mankind realized they could transform it’s surroundings. Furthermore, at that stage the mental maturity achieved made it possible for man to remember how to do it over again in a consistent way.

The book ‘2001 – a space odissey’ by A. C. Clarke, describes this process in the first chapter when a group of hominids bumps into the famous monolith. Their leader (i think his name was moonlight), under the spell of this strangely straight and flat thing takes two pieces of grass and ties them together without knowing or understanding what he is doing. I was pleased to read that I was not alone in that thought.

The concept of a knot keeps on amazing me given their variety and the different purposes they serve according to their properties. These were known to ancient sailors who have elevated the task of knot-making to a practical art form. The mathematical background behind them has served to lay one of today’s most fundamental (and controversial) theories about the composition of matter: string theory. Next time when you make the knot of your necktie think about this tedious, obnoxious little habit was based on something groundbreaking that truly makes us stand out from the rest of the species in the animal kingdom.