# Blog Archives

## Percentage of Molecular Orbital Composition – G09,G16

Canonical Molecular Orbitals are–by construction–delocalized over the various atoms making up a molecule. In some contexts it is important to know how much of any given orbital is made up by a particular atom or group of atoms, and while you could calculate it by hand given the coefficients of each MO in terms of every AO (or basis set function) centered on each atom there is a straightforward way to do it in Gaussian.

If we’re talking about ‘dividing’ a molecular orbital into atomic components, we’re most definitely talking about population analysis calculations, so we’ll resort to the pop keyword and the orbitals option in the standard syntax:

#p M052x/cc-pVDZ pop=orbitals

This will produce the following output right after the Mulliken population analysis section:

Atomic contributions to Alpha molecular orbitals:
Alpha occ 140 OE=-0.314 is Pt1-d=0.23 C38-p=0.16 C31-p=0.16 C36-p=0.16 C33-p=0.15
Alpha occ 141 OE=-0.313 is Pt1-d=0.41
Alpha occ 142 OE=-0.308 is Cl2-p=0.25
Alpha occ 143 OE=-0.302 is Cl2-p=0.72 Pt1-d=0.18
Alpha occ 144 OE=-0.299 is Cl2-p=0.11
Alpha occ 145 OE=-0.298 is C65-p=0.11 C58-p=0.11 C35-p=0.11 C30-p=0.11
Alpha occ 146 OE=-0.293 is C58-p=0.10
Alpha occ 147 OE=-0.291 is C22-p=0.09
Alpha occ 148 OE=-0.273 is Pt1-d=0.18 C11-p=0.12 C7-p=0.11
Alpha occ 149 OE=-0.273 is Pt1-d=0.18
Alpha vir 150 OE=-0.042 is C9-p=0.18 C13-p=0.18
Alpha vir 151 OE=-0.028 is C7-p=0.25 C16-p=0.11 C44-p=0.11
Alpha vir 152 OE=0.017 is Pt1-p=0.10
Alpha vir 153 OE=0.021 is C36-p=0.15 C31-p=0.14 C63-p=0.12 C59-p=0.12 C38-p=0.11 C33-p=0.11
Alpha vir 154 OE=0.023 is C36-p=0.13 C31-p=0.13 C63-p=0.11 C59-p=0.11
Alpha vir 155 OE=0.027 is C65-p=0.11 C58-p=0.10
Alpha vir 156 OE=0.029 is C35-p=0.14 C30-p=0.14 C65-p=0.12 C58-p=0.11
Alpha vir 157 OE=0.032 is C52-p=0.09
Alpha vir 158 OE=0.040 is C50-p=0.14 C22-p=0.13 C45-p=0.12 C17-p=0.11
Alpha vir 159 OE=0.044 is C20-p=0.15 C48-p=0.14 C26-p=0.12 C54-p=0.11


Alpha and Beta densities are listed separately only in unrestricted calculations, otherwise only the first is printed. Each orbital is listed sequentially (occ = occupied; vir = virtual) with their energy value (OE = orbital energy) in atomic units following and then the fraction with which each atom contributes to each MO.

By default only the ten highest occupied orbitals and ten lowest virtual orbitals will be assessed, but the number of MOs to be analyzed can be modified with orbitals=N, if you want to have all orbitals analyzed then use the option AllOrbitals instead of just orbitals. Also, the threshold used for printing the composition is set to 10% but it can be modified with the option ThreshOrbitals=N, for the same compound as before here’s the output lines for HOMO and LUMO (MOs 149, 150) with ThreshOrbitals set to N=1, i.e. 1% as occupation threshold (ThreshOrbitals=1):

Alpha occ 149 OE=-0.273 is Pt1-d=0.18 N4-p=0.08 N6-p=0.08 C20-p=0.06 C13-p=0.06 C48-p=0.06 C9-p=0.06 C24-p=0.05 C52-p=0.05 C16-p=0.04 C44-p=0.04 C8-p=0.03 C15-p=0.03 C17-p=0.03 C45-p=0.02 C46-p=0.02 C18-p=0.02 C26-p=0.02 C54-p=0.02 N5-p=0.01 N3-p=0.01
Alpha vir 150 OE=-0.042 is C9-p=0.18 C13-p=0.18 C44-p=0.08 C16-p=0.08 C15-p=0.06 C8-p=0.06 N6-p=0.04 N4-p=0.04 C52-p=0.04 C24-p=0.04 N5-p=0.03 N3-p=0.03 C46-p=0.03 C18-p=0.03 C48-p=0.02 C20-p=0.02

The fragment=n label in the coordinates can be used as in BSSE Counterpoise calculations and the output will show the orbital composition by fragments with the label "Fr", grouping all contributions to the MO by the AOs centered on the atoms in that fragment.

As always, thanks for reading, sharing, and rating. I hope someone finds this useful.

## Fixing the error: Bad data into FinFrg

I found this error in the calculation of two interacting fragments, both with unpaired electrons. So, two radicals interact at a certain distance and the full system is deemed as a singlet, therefore the unpaired electron on each fragment have opposite spins. The problem came when trying to calculate the Basis Set Superposition Error (BSSE) because in the Counterpoise method you need to assign a charge and multiplicity to each fragment, however it’s not obvious how to assign opposite spins.

The core of the problem is related to the guess construction; normally a Counterpoise calculation would look like the following example:

#p B3LYP/6-31G(d,p) counterpoise=2

-2,1 -1,2 -1,2
C(Fragment=1)        0.00   0.00   0.00
O(Fragment=2)        1.00   1.00   1.00
...

In which the first pair of charge-multiplicity numbers correspond to the whole molecule and the following to those of each fragment in increasing order of N (in this case, N = 2). So for this hypothetical example we have two anions (but could easily be two cations) each with an unpaired electron, yielding a complex of charge = -2 and a singlet multiplicity which implies those two unpaired electrons have opposite spin. But if the guess (the initial trial wavefunction from which the SCF will begin) has a problem understanding this then the title error shows up:

Bad data into FinFrg
Error termination via Lnk1e ...

The solution to this problem is as simple as it may be obscure: Create a convenient guess wavefunction by placing a negative sign to the multiplicity of one of the fragments in the following example. You may then use the guess as the starting point of other calculations since it will be stored in the checkpoint file. By using this negative sign we’re not requesting a negative multiplicity, but a given multiplicity of opposite spin to the other fragment.

#p B3LYP/6-31G(d,p) guess=(only,fragment=2)

-2,1 -1,2 -1,-2
C(Fragment=1)        0.00   0.00   0.00
O(Fragment=2)        1.00   1.00   1.00
...

This way, the second fragment will have the opposite spin (but the same multiplicity) as the first fragment. The only keyword tells gaussian to only calculate the guess wave function and then exit the program. You may then use that guess as the starting point for other calculations such as my failed Counterpoise one.

## Geometry Optimizations for Excited States

Electronic excitations are calculated vertically according to the Frank—Condon principle, this means that the geometry does not change upon the excitation and we merely calculate the energy required to reach the next electronic state. But for some instances, say calculating not only the absorption spectra but also the emission, it is important to know what the geometry minimum of this final state looks like, or if it even exists at all (Figure 1). Optimizing the geometry of a given excited state requires the prior calculation of the vertical excitations whether via a multireference method, quantum Monte Carlo, or the Time Dependent Density Functional Theory, TD-DFT, which due to its lower computational cost is the most widespread method.

Most single-reference treatments, ab initio or density based, yield good agreement with experiments for lower states, but not so much for the higher excitations or process that involve the excitation of two electrons. Of course, an appropriate selection of the method ensures the accuracy of the obtained results, and the more states are considered, the better their description although it becomes more computationally demanding in turn. Figure 1. The vertical excitation does not match the minimum on the excited state

In Gaussian 09 and 16, the argument to the ROOT keyword selects a given excited state to be optimized. In the following example, five excited states are calculated and the optimization is requested upon the second excited state. If no ROOT is specified, then the optimization would be carried out by default on the first excited state (Where L.O.T. stands for Level of Theory).

#p opt TD=(nstates=5,root=2) L.O.T.

Gaussian16 includes now the calculation of analytic second derivatives which allows for the calculation of vibrational frequencies for IR and Raman spectra, as well as transition state optimization and IRC calculations in excited states opening thus an entire avenue for the computation of photochemistry.

If you already computed the excited states and just want to optimize one of them from a previous calculation, you can read the previous results with the following input :

#p opt TD=(Read,Root=N) L.O.T. Density=Current Guess=Read Geom=AllCheck

Common problems. The following error message is commonly observed in excited state calculations whether in TD-DFT, CIS or other methods:

No map to state XX, you need to solve for more vectors in order to follow this state.

This message usually means you need to increase the number of excited states to be calculated for a proper description of the one you’re interested in. Increase the number N for nstates=N in the route section at higher computational cost. A rule of thumb is to request at least 2 more states than the state of interest. This message can also reflect the fact that during the optimization the energy ordering changes between states, and can also mean that the ground state wave function is unstable, i.e., the energy of the excited state becomes lower than that of the ground state, in this case a single determinant approach is unviable and CAS should be used if the size of the molecule allows it. Excited state optimizations are tricky this way, in some cases the optimization may cross from one PES to another making it hard to know if the resulting geometry corresponds to the state of interest or another. Gaussian recommends changing the step size of the optimization from the default 0.3 Bohr radius to 0.1, but obviously this will make the calculation take longer.

Opt=(MaxStep=10)

If the minimum on the excited state potential energy surface (PES) doesn’t exist, then the excited state is not bound; take for example the first excited state of the H2 molecule which doesn’t show a minimum, and therefore the optimized geometry would correspond to both H atoms moving away from each other indefinitely (Figure 2). Nevertheless, a failed optimization doesn’t necessarily means the minimum does not exist and further analysis is required, for instance, checking the gradient is converging to zero while the forces do not. Figure 2. An unbound excited state with no minima ensures the dissociation of the system along the reaction coordinate

## Grimme’s Dispersion DFT-D3 in Gaussian #CompChem

I was just asked if it is possible to perform DFT-D3 calculations in Gaussian and my first answer was to use the following  keyword:

EmpiricalDispersion=GD3

which is available in G16 and G09 only in revision D, apparently.

There are also some overlays that can be used to invoke the use dispersion in various scenarios:

IOp(3/74=x) Exchange and Correlation Potentials

 -77 -76 -60 -59 DSD-PBEP86 (double hybrid, DFT-D3). PW6B95-D3. B2PLYP-D3 (double hybrid, DFT-D3). B97-D (DFT-D3).

IOp(3/76=x) Mixing of HF and DFT.

 -33 PW6B95 and PW6B95-D3 coefficients.

IOp(3/124=x) Empirical dispersion term.

 30 40 50 Force dispersion type 3 (Grimme DFT-D3). Force dispersion type 4 (Grimme DFT-D3(BJ)). Force dispersion type 5 (Grimme D3, PM7 version).

The D3 correction method of Grimme defines the van der Waals energy like: where coefficients are adjusted depending on the geometry of atoms i and j. The damping D3 function for is: where the values of s are adjustable parameters fit for the exchange-correlation functionals used in each calculation.

## Partial Optimizations with Gaussian09

Sometimes you just need to optimize some fragment or moiety of your molecule for a number of reasons -whether because of its size, your current interest, or to skew the progress of a previous optimization- or maybe you want just some kind of atoms to have their positions optimized. I usually optimize hydrogen atoms when working with crystallographic files but that for some reason I want to preserve the rest of the molecule as refined, in order to keep it under a crystalline field of sorts.
Asking Gaussian to optimize some of the atoms in your molecule requires you to make a list albeit the logic behind it is not quite straightforward to me. This list is invoked by the ReadOptimize keyword in the route section and it includes all atoms by default, you can then further tell G09 which atoms are to be included or excluded from the optimization.
So, for example you want to optimize all atoms EXCEPT hydrogens, then your input should bear the ReadOptimize keyword in the route section and then, at the end of the molecule specification, the following line:

atoms notatoms=H

If you wish to selectively add some atoms to the list while excluding others, here’s an example:

atoms=C H S notatoms=5-8

This list adds, and therefore optimizes, all carbon, hydrogen and sulfur atoms, except atoms 5, 6, 7 and 8, should they be any of the previous elements in the C H S list.
The way I selectively optimize hydrogen atoms is by erasing all atoms from the list -using the noatoms instruction- and then selecting which are to be included in the list -with atoms=H-, but I haven’t tried it with only selecting hydrogen atoms from the start, as in atoms=H

noatoms atoms=H

I probably get very confused because I learned to do this with the now obsolete ReadFreeze keyword; now it sometimes may seem to me like I’m using double negatives or something – please do not optimize all atoms except if they are hydrogen atoms. You can include numbers, ranks or symbols in this list as a final line of your input file.

Common errors (by common I mean I’ve got them):

Lets look at the end of an input I just was working with:

>  AtmSel:  Line=”P  0″
>  Maximum list size exceeded in AddBin.
>  Error termination via Lnk1e in…

AtmSel is the routine which reads the atoms list and I was using a pseudopotential on phosphorous atoms, I placed the atoms list at the end of the file but it should be placed right after the coordinates and the connectivity matrix, should there be one, and thus before any external basis set or pseudopotential or any other specification to be read by Gaussian.

As a sort of test you can use the instruction:

%kjob l103
%chk=myfile.chk
...

at the Link0 section (where your checkpoint is defined). This will kill the job after the link 103 is finished, thus you will only get a list of what parameters were frozen and which were active. Then, if things look ok, you can run the job without the %kjob l103 instruction and get it done.

As usual I hope this helps. Thanks for reading except to those who didn’t read it except for the parts they did read.

## Transition State Search (QST2 & QST3) and IRC with Gaussian09

Theoretical evaluation of a reaction mechanism is all about finding the right transition states (TS) but there are no guarantees within the available methods to actually find the one we need. Chemical intuition in the proposal of a mechanism is paramount. Let’s remember that a TS is a critical point on a Potential Energy Surface (PES) that is a minimum in every dimension but one. For a PES with more than two degrees of freedom, a hyper-surface, envisioning the location of a TS is a bit tricky, in the case of a three dimensional PES (two degrees of freedom) the saddle point constitutes the location of the TS as depicted in figure 1 by a section of a revolution hyperboloid. Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

The following procedure considers gas phase calculations. Nevertheless, the use of the SCRF keyword activates the implicit solvent calculation of choice in order to evaluate to some degree the solvent influence on the reaction energetics at different temperatures with the use of the temperature keyword.

The first step consists of a high level optimization of all minimums involved, such as reagents, products and intermediates, with a subsequent frequency analysis that includes no imaginary eigenvalues.

In order to find the structures of the transition states we use in Gaussian the Synchronous Transit-guided Quasi-Newton method  through the keywords QST2 or QST3. In the former case, coordinates for the reagents and products are needed as input; for the latter keyword, coordinates for the TS structure guess is needed also.

QST2)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line

QST3)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title Card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line–
Title card for TS
Q M
Cartesian Coordinates for TS
blank line

NOTE: It is fundamental that the numbering order is kept constant throughout the molecular specifications of all two, or three, input structures. Hence, I recommend to build a set of molecules, save their structure, and then modified the coordinates on the same file to produce the following structure, that way the number for every atom will remain the same for every step.

As I wrote above, there are no guarantees of finding the right TS so many attempts are probably needed. Once we have the optimized structures for all the species involved in our mechanistic proposal we can plot their energies very simply with MS Excel the way we’ve previously described in this blog (reblogged from eutactic.wordpress.com)

Once we’ve succeeded in finding the structure of our TS we may run an Internal Reaction Coordinate (IRC) calculation. This calculation will connect the TS structure to those of the products and the reagents. Initial constant forces are required and these are commonly retrieved from the TS calculation checkpoint file through the RCFC keyword.

%chk=QST3_2p.chk
%nprocshared=8

#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck

Title Card

Q M
blank line

Finally, the IRC path can be visualized with GaussView from the Results menu. A successful IRC will link both structures along a single reaction coordinate proving that both reagents and products are linked by the obtained TS.

Hat tip to Howard Diaz who has become quite skillful in calculating these mechanisms as proven by his recent poster at the XII RMFQT a couple of weeks back. And as usual thanks to everyone who reads, comments, likes, recommends, rates and shares my silly little posts.

## Delta G of solvation in Gaussian09

How to calculate the Delta G of solvation? This is a question that I get a lot in this blog, so it is about time I wrote a (mini)post on it, and at the same time put an end to this posting drought which has lasted for quite a few months due to a lot of pending work with which I’ve had to catch up. Therefore, this is another post in the series of SCRF calculations that are so popular in this blog. For the other posts on this subjects remember to click here and here.

SMD

SMD is the keyword you want to use when performing a Self Consistent Reaction Field (SCRF) calculation with G09. This keyword was only made available in this last version of the program and it corresponds to Truhlar’s and coworkers solvation model which is recommended by Gaussian itself as the preferred model to calculate Delta G of solvation. The syntax used is the standard way used in any other Gaussian input files as follows:

# 'route section keywords' SCRF=SMD

Separately, we must either perform a gas phase calculation or use the DoVacuum keyword within the same SCRF input, and then take the energy difference between gas phase and solvated models.

# 'route section keywords' SCRF=(SMD,DoVacuum)

No solvation or cavity model should be defined since, by definition, SMD will use the IEFPCM model which is a synonym for PCM.

As opposed to the previous versions of Gaussian, the output energy already contains all corrections, this is why we must take the difference between both values (remember to calculate them both at the same level of theory if calculated separately!). Nevertheless, when using the SMD keyword we get a separate line, just below the energy, stating the SMD-CDS non electrostatic value in kCal/mol.

The radii were also defined in the original paper by Truhlar; I’m not sure if using the keyword RADII with any of its options yields a different result or if it even ends in an error. Its worth the try!

Some calculation variations are not available when using SMD, such as Dis (calculation of the solute-solvent dispersion interaction energy), Rep (solute-solvent repulsion interaction energy) and Cav (inclusion of the solute cavitation energy in the total energy). I guess the reason for this might be that the SMD model is highly parametrized.

Have you found any issue with any item listed above? Pleases share your thoughts in the comments section below. As usual I hope this post was useful and that you all rate it, like it and comment.

References

A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Universal solvation model based on solute electron density and a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions,” J. Phys. Chem. B, 113 (2009) 6378-96.