Eduticles » Molecular Modeling Basics

« Expand/Collapse

Molecular Modeling Basics (4 unread)

Skip to page: 1 2 3

March 03, 2010
2:58
Internal energy and molecular motion
» ‎ Molecular Modeling Basics

This screencast shows a Molecular Workbench model I made to illustrate the connection between internal energy and molecular motion (the model is included at the bottom of the post).

For a monatomic gas such as He or Ar the part of the internal energy (U) that depends on temperature (T) is
and has only contributions from translational motion. Kinetic theory tells us that
which, with a little algebra, gives us a direct connection between the internal energy and the average kinetic energy of the atoms in the gas:
You can see in the first container that while the internal energy reflects an average speed, the individual atoms can have very different speeds, which are constantly changing due to collisions. This can be illustrated by coloring the atoms according to the kinetic energies (red means more kinetic energy).

Increasing the temperature increases the internal energy and this is reflected in an increase in the average speed, as you can see by clicking on the "Increase T" button.

However, it is important to remember that the internal energy is a reflection of the kinetic energy, which is also a reflection of the mass. Thus, a heavier atoms will have the same internal energy and kinetic energy at the same temperature, but the atoms will move slower on average (the "Increase M" button).

Also, the internal energy can be increase by adding more particles, but this does not change the average speed of the molecules at the same temperature (the "Increase n" button).

The last ("Diatomic") button changes the atoms to diatomic molecules without changing the mass, i.e. the mass of the atoms in the second container are half that of the atoms in the first container. Because the mass is unchanged the average speed of the molecules (more precisely, of their center of mass) is the same as for the atomic gas at the same temperature, as is the translational internal energy (³/2nRT ).

However, notice that the molecules are also rotating, meaning they have additional (internal) energy compared to the atomic gas, which amounts to nRT (³/2nRT for non linear molecules):
The internal energy also has a vibrational contribution (here the frequency is converted to a wavenumber, i.e. units of cm^-1):
though it is hard to make out the vibrational motion in the simulation. However, for most diatomic molecules this contribution is negligible at temperatures relevant to most chemists.

You can play around with the model yourself by clicking on the picture below, or on this web page, or you can download the model if you have Molecular Workbench installed on your computer. Enjoy!
February 20, 2010
14:32
Ionic and metallic bonding
» ‎ Molecular Modeling Basics

Figure 4.7. 0.0005 au isodensity surface with superimposed electrostatic potential of (a) H–Li, (b) Li, (c) Li2. The maximum potential value is 0.05 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.
Note: the figure is different from the book. See details at the end of the post.

Figure 4.7 highlights two other kinds of bonds. Note the decrease in the size of the Li atom (and the large positive charge) compared to the Li atom in Li–H. LiH is best thought of as Li⁺H^- (i.e., the bonding electron pair belongs to H^- rather than being shared), which is an example of ionic bonding.

Here is an interactive version of Figure 4.7a where the atomic densities of H and Li are superimposed on the density of H-Li.

Click on the picture for an interactive version

Li2 is more polar than H2: more positive at the ends and more negative in the middle, indicating the a larger rearrangement of electron density upon binding compared to H2. Here is an interactive version (the corresponding Jmol script can be found here; see this post on making plots like this), where I have superimposed the electron densities of the Li atoms and Li2:
Click on the picture for an interactive version

Looking deeper into the density (by using a larger, 0.01 au, isodensity value) reveals a density rearrangement that is quite different from H2:
Click on the picture for an interactive version

Here you can see strong localization of electron density between the nuclei, indicating that Li2 is best thought of as Li⁺:Li⁺, where ":" indicates an electron pair. This is an example of metallic bonding.

Both ionic and metallic bonding are distinctly different from covalent bonding, in that only covalent bonding leads to the formation of distinct molecules, while ionic and metallic bonding leads to formation of crystals.

To demonstrate this, find the lowest energy structures of H4, Li4, and Li2H2, using, for example, B3LYP/6-31G(d). Go ahead, I'll wait.

How and why is the figure different from the book?
When prepring this post I discovered that I had shown an H atom instead of a Li atom in Figure 4.7b. When I made the plot of the 0.002 au isodensity surface with superimposed electrostatic potential of Li atom, I discovered that is was very blue, i.e. quite positive (go ahead, try it for yourself). I only got a neutral-looking atom when I went down to 0.0005 aus, which is what I have used for the figures (unless otherwise noted).

The positive electrostatic potential on the 0.002 au isodensity surface of Li indicates that a significant amount of density is found outside this surface, which is clearly not the case of the H atom, as shown in a previous post. This, in turn, indicates that the Li electron density decreases more slowly with distance compared to H, which makes sense since Li is much less electronegative than H. Thus, the 0.002 au isodensity surface is a reflection of the van der Waals surface of organic molecules (containing electronegative elements) but not alkali metals.
February 18, 2010
14:14
Covalent bonding
» ‎ Molecular Modeling Basics

Figure 4.5. 0.002 au isodensity surface with superimposed electrostatic potential of (a) H2 and (b) H atom. The maximum potential value is 0.05 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

If we look at the electron density of H2 (Figure 4.5a), we can clearly see that at this separation the electron densities of the two H atoms have fused indicating electron sharing, a hallmark of covalent bonding. Here is an interactive version (the corresponding Jmol script can be found here; see this post on making plots like this), where I have superimposed the electron densities of the H atoms and H2:

Click on the picture for an interactive version

What you can't see in this picture is that the electron density has rearranged significantly between the nuclei, which is the source of the bond strength. To see this we need to look at a larger isodensity value (here 0.075 au):

Click on the picture for an interactive version
February 16, 2010
1:58
One, two, three, MD
» ‎ Molecular Modeling Basics

Some exciting developments over at the Molecular Workbench (MW) blog run by MW author Charles Xie. I have several blog posts on MW, but it was necessary to install MW to check it out for yourself. Now it is possible to embed a MW applet in web pages (and blog posts!) like the one here (just push the play button!):

I think this is a big step forward for MW. While it is easy to download and install MW, it still removed MW a few clicks from the user and made it "appear to be yet another kind of annoying pop-up" and Charles notes.

It's very easy to do this. The screencast below shows how I made the simulation above in MW. Note that it literally takes one minute (and 5 seconds).

When you hit save you get two files: md.cml and md$0.mml. I transferred these to my web server where I had also put the MW applet (mwapplet.jar).

The html code is

<applet code="org.concord.modeler.MwApplet" archive="mwapplet.jar" height="300" width="100%"><param name="script" value="page:0:import md.cml"></applet>

To include it in a blog, where mwapplet.jar is not installed, add the server address in front of mwapplet.jar and md.cml, e.g. [myserver.edu]
February 07, 2010
8:57
I'm positive
» ‎ Molecular Modeling Basics

Figure 4.2. 0.002 au isodensity surface with superimposed electrostatic potential of (a) Li⁺, (b) Na⁺, and (c) K⁺ ion. The maximum potential value is 0.8 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

Here is an example of how computational chemistry can be used to enhance teaching at the general chemistry level.

Why does the ionization energy decrease on from Li to Na to K? That's the same as asking why the electron affinities decrease on going from Li⁺ to Na⁺ to K⁺.

Figure 4.2 show the ions colored by how positive they are at the surface (i.e. the electrostatic potential superimposed on the 0.002 isodensity surface). The darker the
color the more positive the ion, and it is clear that the Li+ ion is “more positive” than Na⁺, which is more positive than K⁺. Thus, more energy should be released when adding an electron to Li⁺ compared to Na⁺, and hence more energy is needed to remove an electron from Li
compared to Na (and similarly for K).

The reason why Li⁺ is “more positive” than Na⁺ or, more accurately, why the potential on the 0.002 au isodensity surface is more positive for Li⁺ than for Na⁺, is that the former is a smaller ion than the latter, so the surface is closer to the +1 charge at the center of the ion.

Figure 4.4. 0.002 au isodensity surface with superimposed electrostatic potential of (a) Li⁺, (b) Ne⁺, and (c) Na⁺ ion. The maximum potential value is 0.8 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

This rationalization is, of course, only qualitative and is is not predictive of the ionization energies among different groups. For example, based on Figure 4.4 one would expect that Ne would have roughly the same ionization potential as Na, which is not true at all.

When making these figures it is very important to get the relative sizes of the ions correct, but this can be difficult since each image can be zoomed to an arbitrary size. The screencast below shows how to control this in MacMolPlt using the Manual Windows Parameter window.

The same point can also be made with a "quick-and-dirty" electrostatic potential, as shown in this interactive figure. Here the electrostatic potential is due to a plus one charge centered at the atom (rather than the nuclear charge and the electron density) and the surfaces are the spheres defined by empirical ionic radii. The Jmol script can be found here, and the mol2 files here, here, and here.

Click on the picture for an interactive version
January 31, 2010
2:27
Tesserae, sera
» ‎ Molecular Modeling Basics

Figure 3.14. A solute–solvent boundary surface described by a set of interlocking atom-centered spheres. The surface is discretized by using triangular tesserae. Drawing graciously provided by Benedetta Mennucci.
From Molecular Modeling Basics CRC Press, May 2010.

Due to the computational expense of explicit solvation, implicit solvation is a popular alternative when including solvent effects in quantum chemical calculations. The polarized continuum model (PCM) and the closely related COSMO method is one of the most used implicit solvation models for QM calculations. Here the Poisson-Boltzmann equation is solved numerically by diving the solute-solvent surface into small pieces called tesserae (Figure 3.14).

Unfortunately, I have not been able to make an interactive version of this figure.
January 27, 2010
15:08
It's all about boundaries
» ‎ Molecular Modeling Basics

Figure 3.8. (a) Adenine (C5H5N5) in a liquid drop of 246 water molecules. (b) Adenine in a periodic box of 511 water molecules.
From Molecular Modeling Basics CRC Press, May 2010.

Conceptually, the simplest way of simulating a molecule in solution is to place it in the middle of a roughly spherical ball of water molecules (Figure 3.8a) and perform an molecular dynamics simulation.

One problem is this approach is that the drop would eventually evaporate if the simulation is run long enough. Another problem with the liquid drop model is that the water molecules at the surface of the drop do not behave like water molecules in liquid water.

Therefore, most explicit solvation simulations use periodic boundary conditions (Figure 3.9).

Figure 3.9. Sketch of periodic boundary in two dimensions: (a) The position of the particles in the central box are copied and placed in neighboring boxes. Figure 3.8b shows a cube from a real simulation. (b) When a molecule tries to leave the box during an MD simulation, it reappears at the opposite end of the box, so the number of particles in the central box stays constant.
From Molecular Modeling Basics CRC Press, May 2010.

You can find interactive versions of Figures 3.8a and 3.9b here (I am grateful to Kestutis Aidas for providing the coordinates).

Click on the picture for an interactive version

Click on the picture for an interactive version

And you can find an animated version of Figure 3.9 here (an example of where a movie really is worth 10,000 words).

The animation was made with Molecular Workbench (MW). You can download the MW file here (after you gave installed MW).
January 23, 2010
6:28
I have my moments
» ‎ Molecular Modeling Basics
Figure 3.7. Contour plot of the RHF/6-31G(d) electrostatic potential and 0.002 au isodensity surface of (a) CH3COO^-, (b) HF, and (c) F2. The maximum/minimum contour values are, respectively, 0.5/0.025; 0.1/0.005; and 0.005/0.00025 au respectively. Blue corresponds to a negative potential. In each case the outer-most contour looks like the corresponding contour in the electrostatic potential due to a charge, dipole, and quadrupole, respectively.
From Molecular Modeling Basics CRC Press, May 2010.

This figure makes 3 points:

1. It shows what the electrostatic potential of a charge, dipole, and quadrupole looks like.

2. It show the relative strengths of the electrostatic potentials due to a charge, dipole, and quadrupole.

3. It shows that the electrostatic potential of a charge, dipole, and quadrupole deviates significantly from the actual electrostatic potential near the molecular surface.

Here is a screencast showing how I made Figure 3.7a. (If I were to do it over I would have chosen red for negative and blue for positive... oh, well.)

Here is an interactive version of the figure. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.7. The RHF/6-31G(d) electrostatic potentials of acetate, HF and F2.
Click on the picture for an interactive version

See these two posts (here and here) on how to make plots like that with Jmol. From a Jmol perspective the only new thing is that I show two isosurfaces simulateneously. This is done by loading the file twice to create two "frames" that Jmol can display simultaneously. You can find the script file with all the commands here, but the general syntax is:
```
load files file1.xyz file2.xyz
frame 1.1; isosurface surf1 plane {0 0 0 0} contour 40 color range -0.05 0.05 "potential.cube.gz"
frame 2.1; isosurface surf2 0.002 "density.cube.gz" 
```
Even though you want a 2D contour plot of the potential, it is necessary to make a 3D cube file. Because I use unusually low cutoffs to show the outer contours, I had to trick MacMolPlt into making a bigger grid. I show how on the screencast below.
January 02, 2010

Permalink for 'Common error messages in GAMESS: Failure to locate stationary point'

11:45

Common error messages in GAMESS: Failure to locate stationary point

» ‎ Molecular Modeling Basics

       ***** FAILURE TO LOCATE STATIONARY POINT, TOO MANY STEPS TAKEN *****
  UPDATED HESSIAN, GEOMETRY, AND VECTORS WILL BE PUNCHED FOR RESTART
**** THE GEOMETRY SEARCH IS NOT CONVERGED! ****

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize icharg=1 $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

Actually, there is no problem with the input file as such. A geometry optimization is an iterative process and if the gradient it not below the convergence criteria within 20 steps, GAMESS will stop and print out the message shown above.

The solution is simply to take the last set of coordinates and run the optimization again, as I show in the screencast below.

As I've mentioned in a previous post I think the default criterion for geometry convergence (0.0001) is too strict, and the default number of steps (20) is too small. So I usually use 0.0005 and 50, respectively.

 $statpt nstep=50 opttol=0.0005 $end

December 30, 2009
14:21
Quick and dirty electrostatic potential maps
» ‎ Molecular Modeling Basics
In a previous post I showed how to compute an electrostatic potential map superimposed on the 0.002 isodensity surface of a molecule based on data computed using quantum chemical methods such as RHF/6-31G(d).

Previous posts (such as this one) have demonstrated that the van der Waals surface is a reasonable substitute for the 0.002 isodensity surface. Similarly, the electrostatic potential due to atomic charges can be a reasonable substitute for the electrostatic potential due to the electronic density.

The screencast above shows how to make such electrostatic potential maps using Avogadro and Jmol.

Avogadro uses an empirical method to determine the atomic charges (an integral part of the MMFF force field). It is possible to change the surface using the so-called "Iso Value" but I could not find any documentation on how that actually works. An Iso Value of 0 seems to correspond to the van der Waals sphere surface. It is currently not possible to alter the color range as far as I can see.

Jmol is not able to determine the charges, but the information can be transferred from Avogadro by saving a mol2 file. The electrostatic potential option in the Jmol menu corresponds to the following set of commands (as far as I can determine):
```
isosurface solvent color range -0.05 0.05 map mep
color isosurface translucent 0.5
```
and this set of commands can thus be used to control the color range and the nature of the surface. The default surface is a solvent accessible surface, which is slightly larger than the van der Waals surface.

A static and interactive version of the Avogadro and Jmol electrostatic potential map, respectively, can be found here.

Electrostatic potential map made with Avogadro
Click on the picture for an interactive version made with Jmol
December 29, 2009
8:37
Steric strain vs electrostatic attraction
» ‎ Molecular Modeling Basics

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH3(H)C=C(H)CH3 and. In both cases, the maximum potential value is 0.05 au. (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

This figure shows how the electrostatic potential superimposed on the 0.002 au isodensity surface can be used to rationalize why cis-HO(H)C=C(H)OH is more stable than trans-HO(H)C=C(H)OH, while the opposite is true for CH3(H)C=C(H)CH3.

The figure clearly shows the difference in polarity between hydroxyl and methyl groups. For the OH substitutent case the positive (O)H atom is close to the negative O(H) atom in the cis isomer. The CH3 group is non-polar and larger and the overlapping density indicates steric strain.

The figure is made with MacMolPlt as described in a previous post. Below is an interactive version made with Jmol (as described in a previous post). You'll notice that the color scheme is somewhat different, because Jmol maps the electrostatic potential value to a color in a different way than MacMolPlt. However, the general conclusion drawn from both programs is clearly the same.

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH3(H)C=C(H)CH3 and.
Click on the picture for an interactive version
December 27, 2009
2:58
Electrostatic potential maps reloaded
» ‎ Molecular Modeling Basics
Here is an interactive version of a figure I described in a previous post. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.4. The RHF/6-31G(d) electrostatic potential of water.
Click on the picture for an interactive version

The Jmol animation loads a file (fig3-2-2.xyz), which you can access by right-clicking on the Jmol animation once you have loaded it (as described in a previous post). This file also contains all the necessary commands.

Figure 3.4.c is the most common depiction of electrostatic potential maps and the Jmol general syntax for the command for this is
```
isosurface 0.002 "density.cube.gz" color range -0.05 0.05 "potential.cube.gz"
```
The screencast below shows how I made the cube files that contain the electron density and electrostatic potential information using MacMolPlt. The first part is identical to the screencast in a previous post on obtaining the electron density cube file.

The program I used to convert the MacMolPlt file to a cube was written by Jonathan Gutow and can be download here.

I start by loading the RHF/6-31G(d) optimized geometry I computed for a previous post, and reorienting. The orientation makes it easier to define the plotting plane for the contour plots.
December 26, 2009
3:25
Wishlist
» ‎ Molecular Modeling Basics

Christmas is only 363 days away. It's never to early to start working on a new wish list.

[source]
[source]
[source]

[source]
December 21, 2009
2:35
Electrostatic potential maps
» ‎ Molecular Modeling Basics

Figure 3.4. (a) Contour plot of the electrostatic potential of H2O. The maximum and minimum contour value is 0.5 Hartrees/electronic charge (au). (b) The corresponding 0.05 au isopotential surface. (c) The electrostatic potential displayed on the 0.002 au isodensity surface of water. The maximum (darkest blue) value corresponds to 0.05 au.
From Molecular Modeling Basics CRC Press, May 2010.

Here is a screencast of how I made the figure:

The files I used were created in a previous post. The various values specified above are determined using trial and error, i.e. I kept fiddling with it until it looked good to me. When using plots like this it is important to specify these values, because using different values can lead to very different looking plots.

Also, different programs use different color scales and color intensity to indicate positive and negative charge, so it is rarely possible to directly compare electrostatic potential maps from different programs.

See this post about using screen capture to get a file with the graphic. This was done for each plot, and the combined in a word processor.

The interactive version of this figure is the subject of a future post.
December 19, 2009
4:00
Electron density doesn't always tell the whole story
» ‎ Molecular Modeling Basics

Figure 3.3. (a) RHF/6-31G(d) 0.002 au isodensity surface and (b) van der Waals surface of cis-HO(H)C=C(H)OH (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

Looking at the electron density or the van der Waals surfaces one would expect steric strain in cis-HO(H)C=C(H)OH. (See this post on how to make such pictures and this post on how to make an interactive version such as the one shown below.)

However, contrary to CH3(H)C=C(H)CH3, the cis isomer is more stable than the trans isomer. This is of course due to the electrostatic attraction of the O and H atom in the hydrogen bond. Here the electron density tells only part of the story because this attraction involves the nuclei as well. Here an electrostatic potential map is more useful. How to make such plots will be the subject of future posts.

Figure 3.3. The RHF/6-31G(d) electron density of cis-HO(H)C=C(H)OH.
Click on the picture for an interactive version
December 11, 2009

Permalink for 'Common error messages in GAMESS: charge and multiplicity'

7:05

Common error messages in GAMESS: charge and multiplicity

» ‎ Molecular Modeling Basics

  *** CHECK YOUR INPUT CHARGE AND MULTIPLICITY ***
THERE ARE    33 ELECTRONS, WITH CHARGE ICHARG=  0
BUT YOU SELECTED MULTIPLICITY MULT=  1

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

As GAMESS tells you, it is not possible to have singlet state (i.e. a multiplicity of 1) if you have an odd number of electrons. You need to tell GAMESS whether you want to change the multiplicity (for example, $contrl mult=2 scftyp=rohf $end) or add ($contrl icharg=-1 $end) or remove ($contrl icharg=1 $end).

When looking at the molecule it is pretty clear that you want to remove an electron, since ammonium has a positive charge. So the first line in the input should be changed to

 $contrl runtyp=optimize icharg=1 $end

Click on the picture for an interactive version

Notice that you don't have to specify that the positive charge is on ammonium. The quantum mechanics will take care of that.

December 09, 2009
6:58
Electron density and steric strain
» ‎ Molecular Modeling Basics

Figure 3.2. (a) RHF/6-31G(d) 0.002 au isodensity surface and (b) van der Waals surface of cis-CH3(H)C=C(H)CH3 (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

This figure shows how the electron density and van der Waals surfaces can be used to visualize steric strain. See this post on how to make such pictures and this post on how to make an interactive version such as the one shown below.

Figure 3.2. The RHF/6-31G(d) electron density of cis-CH3(H)C=C(H)CH3.
Click on the picture for an interactive version
December 08, 2009

Permalink for 'Common error messages in GAMESS: error reading IDUM'

14:13

Common error messages in GAMESS: error reading IDUM

» ‎ Molecular Modeling Basics

 **** ERROR READING VARIABLE IDUM     CHECK COLUMN  1
H     1.0    -0.08386     2.76975     0.70945                                
....V....1....V....2....V....3....V....4....V....5....V....6....V....7....V....8

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize $end
$basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $basis, and the solution is to add a space in front of it, just like for $contrl and $basis.

Why does GAMESS complain about a line in $data when the problem is with $basis? The reason is that it is possible to specify a basis set for each atom in $data, so GAMESS tries to read this information if no $basis group is found.

The error message would clearly benefit from the line: "Did you mean to specify a $basis group?"

December 06, 2009
4:00
The future of scientific publishing is here today
» ‎ Molecular Modeling Basics

I recently came across this paper on embedding interactive 3D graphics in pdf files. Unfortunately, I couldn't animate the 3D figures, and I wrote the author (Vlad Vasilyev) who kindly wrote back that:

(1) the publisher (Springer) is a bit behind the times so the interactive figures are found in the supplemental material, but

(2) the ACS is more on the ball, as can be seen in this paper (this is the one I use in the screencast), and

(3) he has made a bunch of tutorials on how to create such diagrams on this page.

This page also contains some pdf files with 3D graphics you can play around with if you don't have access to the papers. Be sure to use Adobe Acrobat Reader 9.

Note that it is also possible to include animations such as vibrational motion (NB: 16 MB file)

Very cool! I haven't tried this myself, but you can be sure I will. Unfortunately, it seems that you need a piece of commercial software (Adobe Acrobat 9 Pro Extended) to do this.

On a related note, check you the interactive figure of interactive figures created by fellow blogger Henry Rzepa (please give the page some time to load).
November 29, 2009
8:26
Electron density reloaded
» ‎ Molecular Modeling Basics
Here is an interactive version of a figure I described in a previous post. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.1. The RHF/6-31G(d) electron density of water.
Click on the picture for an interactive version

The Jmol animation loads this file (h2oprinc.xyz), which looks like this
```
3
jmolscript: script "http://propka.ki.ku.dk/~jhjensen/h2odensity.spt"
O -0.0000 0.0643 0.0000
H -0.7541 -0.5091 0.0000
H 0.7541 -0.5091 0.0000
```
and which, in turn, loads a script (h2odensity.spt), which looks like this
```
isosurface planex plane {0 0 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
delay 3
spin y 20
delay 10
spin off
isosurface planey plane {1 0 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
spin y 20
delay 10
spin off
isosurface planez plane {0 1 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
spin x 40
delay 10
spin off
isosurface threed 0.002 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
color isosurface red ; color isosurface translucent 0.15
spin y 20
delay 10
spin off
color isosurface red ; color isosurface translucent 0.5
select all; spacefill 100 %babel
spin y 20
delay 10
spin off
```
The screencast below shows how I made the cube file that contains the electron density information using MacMolPlt. The program I used to convert the MacMolPlt file to a cube was written by Jonathan Gutow and can be download here. In the screencast I mistakenly named the cube file h2odensity.cube.gz, but that's easy to change.

I start by loading the RHF/6-31G(d) optimized geometry I computed for a previous post, and reorienting. The orientation makes it easier to define the plotting plane for the contour plots.
November 25, 2009

Permalink for 'Common error messages in GAMESS: No $data'

1:36

Common error messages in GAMESS: No $data

» ‎ Molecular Modeling Basics

**** ERROR, NO  $DATA   GROUP WAS FOUND
EXECUTION OF GAMESS TERMINATED -ABNORMALLY-

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
$data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $data, and the solution is to add a space in front of it, just like for $contrl and $basis.

November 21, 2009
6:00
Electron density
» ‎ Molecular Modeling Basics

Figure 3.1. (a) A contour plot of the density of H2O computed using RHF/6-31G(d). The maximum and minimum contour values are 0.05 and 0.002 aus. (b) The corresponding 0.002 au isodensity surface. (c) The surface corresponding defined by atomic spheres with van der Waals radii.
From Molecular Modeling Basics CRC Press, May 2010.

Here's a screencast about how I made the figure:

When making the contour plot (Figure 3.1.a) I pick 0.05 au as the maximum value (this will be the contour line closest to the nuclei) and 25 contour lines. This means that the spacing between the contour lines and thus the outer contour line will be 0.05/25 = 0.002 au.

Another, more common, representation of the density is a 3D version of one of the contour values: the isodensity surface (Figure 3.1.b). A common choice is 0.002 au, since that corresponds roughly to experimental estimates of molecular size, such as the van der Waals surface (Figure 3.1.c).

Unfortunately, MacMolPlt doesn't have a van der Waals display style, so I have to use Avogadro. This means I have to re-size (by eyeball) this part of the figure to make it the same size as the density plots. See this post about using screen capture to get a file with the graphic.

The interactive version of this figure is the subject of a future post.
5:18
The book and color figures
» ‎ Molecular Modeling Basics

A big part of the motivation for this blog came from writing a book called Molecular Modeling Basics that will be published in May, 2010 by CRC Press. While writing the applications sections it was frustrating to turn the beautifully colored figures into black-and-white versions in order to keep the cost of the book reasonable. But it was also apparent that even colored figures in a book would be a somewhat poor substitute for the interactive versions they are based on. Especially, when turning them around to find just the right orientation for the figure. Wouldn't it be much better to have the reader decide for him/herself?

This is all a long winded way of explaining why there'll be a lot of posts with (color) figures that look like they came out of a book (they'll have figure captions below them). You can click on them for a bigger version. In many of the posts there'll also be a screencast showing how I made them, and an interactive Jmol version. They'll all be labelled "color figures from the book" so they should be easy find.

I always use a screen capture program, rather than saving a graphics file (see the screencast below for an example). I use a free Mac program called Snap and Drag, but I suspect there are plenty of other options for Windows and Linux.
November 11, 2009
14:00
A sense of scale
» ‎ Molecular Modeling Basics

In chapter 1 of volume 1 of the legendary Feynman Lectures on Physics, Feynman starts by imagining zooming in on a drop of water, past amoeba and so forth, until one can see the water molecules. Starting this way is genius. The tiny length scales and the associated invisibility of atoms and molecules is the single largest barrier to developing a chemical intuition - a barrier that molecular animation can help overcome.

I am often thought of bringing Feynman's imaginary magnification to life by animation, so I as very happy when Nathan Baker brought this site to my attention. The above screencast shows it in action, but the real fun is interacting with it yourself. It's a brilliant piece of interactive animation.

The site brought to mind the granddaddy of them all, Powers of 10, which some kind soul put up on youtube.
November 09, 2009

Permalink for 'The trouble with transition metals II: SCF convergence'

15:03

The trouble with transition metals II: SCF convergence

» ‎ Molecular Modeling Basics

In a previous post I showed how to build two structural models containing transition metals: a small model of a zinc finger and a ferrocene. This post is about computing the energy using GAMESS.

The first challenge in computing the energy is choosing the correct charge and multiplicity. Most organic molecules are singlets (i.e. all orbitals are doubly occupied) and have an easily identifiable charge. This is often not the case with molecules containing transition metals. Many transition metal atoms have unpaired d electrons and can have several different charges (oxidation states). Also, the charge of the ligands is not always obvious, though they almost always are singlets.

The zinc finger model is relatively easy: zinc is almost always Zn²⁺ and has 5 doubly occupied d orbitals, so that's a pretty safe bet. The zinc finger model contains two neutral groups (the imidazole rings) and two negatively charged groups (CH3-S^-). So the charge of the entire system is zero and the multiplicity is 1.

The ferrocene model is a bit more complicated: iron is commonly found in both the +2 and +3 oxidation state. Before one starts the calculation it is important to find out which one is appropriate for ferrocene. Google says +2 (meaning 6 d electrons) and all doubly occupied orbitals (see the orbitals at the bottom of the page), so a singlet. Furthermore, it is also important to know that each ring has a -1 charge, so the total charge of the complex is 0 and the mutiplicity is 1.

Now, I was all set to talk about SCF convergence problems, i.e. that the zinc finger would converge fine but the ferrocene would give problems. But when I ran a RHF/6-31G(d) energy calculation

$contrl runtyp=energy $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$scf dirscf=.t. $end

both runs converged fine!

"Luckily" when using a smaller basis set (RHF/STO-3G) leads to convergence problems - for both. Here is the relevant part of the SCF output for zinc finger.

                                                                                                                                                                     NONZERO     BLOCKS
ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE     ORB. GRAD      INTEGRALS    SKIPPED
 1  0  0    -3064.6703808712 -3064.6703808712   1.733010527   0.000000000        9751402     457807
 2  1  0    -3059.0602240553     5.6101568159  16.724847696   0.794728909        9737106     471195
 3  2  0    -3053.0258026844     6.0344213709  16.846293078   0.656307634        9778274     459644
 4  3  0    -2935.9005712044   117.1252314800 115.104028436   1.822663152        9792060     446482
 5  4  0    -2923.8944526114    12.0061185930 115.740577390   0.620923864        9792989     441267
 6  5  0    -2768.2361993103   155.6582533011 115.747362137   5.442595435        9783876     443565
 7  6  0    -2921.1390148151  -152.9028155048 115.772970667   0.810505145        9782321     444765
 8  7  0    -2775.0326048138   146.1064100013 115.770983819   4.741230605        9783214     444582
 9  0  0    -2920.9171138162  -145.8845090024  60.570511833   0.796143103        9783475     444487
10  1  0    -3003.9501726068   -83.0330587906  54.636432025   0.749401761        9790866     444493
11  2  0    -2862.1853212282   141.7648513786 115.095938299   1.824897836        9796139     439522
12  3  0    -2921.4368542174   -59.2515329892 115.495602218   0.722894426        9783664     443903
13  4  0    -2848.2964423996    73.1404118178 115.538938253   2.394847457        9782038     444795
14  5  0    -2905.4155965583   -57.1191541586 115.529755578   0.678236838        9781255     445351
15  0  0    -2859.1068450197    46.3087515386  14.660262004   1.964358011        9783437     444643
16  1  0    -3019.7713654862  -160.6645204665  14.032213984   0.801579670        9794434     443563
17  2  0    -2934.9445251088    84.8268403773  99.891887745   1.748565491        9794346     445610
18  3  0    -2905.6599083961    29.2846167127 100.392269592   1.167341691        9789953     445157
19  4  0    -2864.9366576173    40.7232507788 104.001625891   1.797814244        9781805     446591
20  5  0    -2855.1081148594     9.8285427580 103.862237046   1.447532031        9781260     446026
21  6  0    -2865.3401581051   -10.2320432457  82.200569407   1.143844678        9780463     446052
22  0  0    -2857.7261418857     7.6140162194  78.128990865   1.766124007        9781544     445447
23  1  0    -3036.7926259046  -179.0664840190   5.185277571   0.905312779        9796083     442344
24  2  0    -2977.1610659805    59.6315599241   5.419478057   1.193389134        9796601     444172
25  3  0    -2859.5939849458   117.5670810347 114.799971576   1.724872824        9789000     444821
26  4  0    -2913.7281172708   -54.1341323250 115.553208606   1.094842493        9785096     443663
27  0  0    -2845.9860251842    67.7420920865  16.890315601   1.995942532        9783715     444171
28  1  0    -2989.2091196251  -143.2230944409  16.260193167   1.244487455        9796844     439846
29  2  0    -2911.6240721064    77.5850475187 114.870828357   1.711669748        9798193     439590
30  3  0    -2897.8498052725    13.7742668339 115.443658127   0.704588469        9789911     442464

SCF IS UNCONVERGED, TOO MANY ITERATIONS
  TIME TO FORM FOCK OPERATORS=      84.7 SECONDS (       2.8 SEC/ITER)
  FOCK TIME ON FIRST ITERATION=       2.6, LAST ITERATION=       3.0
  TIME TO SOLVE SCF EQUATIONS=       0.4 SECONDS (       0.0 SEC/ITER)

FINAL RHF ENERGY IS        0.0000000000 AFTER  30 ITERATIONS

Clearly, the this SCF is not converging and giving it more iterations is not going to make the problem go away. Instead, different algorithms for SCF convergence are needed and

$scf dirscf=.t. diis=.t. damp=.t. $end

leads to convergence in for both ferrocene and zinc finger.

DIIS stands to Direct Inversion of the Iterative Subspace, which is an extrapolation procedure, while DAMP "damps" any oscillations between in the SCF.

There are a few other tricks to SCF convergence but they are clearly not needed here. If DIIS and DAMP doesn't solve the problem for you, leave a message describing your molecule and I'll see what I can do.

Skip to page: 1 2 3

← Phosphorus carbide Indiana University Cheminformatics News →

Eduticles » Molecular Modeling Basics

Feeds

Molecular Modeling Basics (4 unread)